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ABSTRACT 


The  analysis  of  multiple  correlated  two-dimensional  (2-D)  random  signals  or  mul¬ 
tichannel  2-D  signals  is  described.  The  emphasis  is  on  estimation  (linear  prediction) 
and  modeling  of  the  2-D  random  signals.  Applications  to  spectrum  analysis  and  image 
processing  are  considered. 
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I.  Introduction 


This  report  deals  with  the  analysis  of  multichannel  two-dimensional  (2-D)  random 
signals.  The  signals  are  two-dimensional  in  that  the  signal  is  a  function  of  two  inde¬ 
pendent  variables  {i%i ,  n2).  We  say  the  signals  are  multichannel  since  there  are  possibly 
several  correlated  2-D  signals,  perhaps  received  on  different  communication  channels 
that  are  processed  and  analyzed  together. 

One  example  of  a  multichannel  2-D  signal  is  the  set  of  images  received  from  a 
satellite  multispectral  scanner.  The  satellite  forms  several  images  of  the  ground  in  pixel 
registration  but  corresponding  to  different  frequency  bands  in  the  visible  and  infrared. 
One  can  have  four,  seven,  or  even  more  such  channels  comprising  the  multichannel 
data. 

A  more  common  example  of  a  multichannel  2-D  signal  is  a  color  image.  The  signal 
is  inherently  two-dimensional  but  consists  of  three  registered  components  representing 
red,  green  and  blue  intensities  or  other  tristimulus  values  such  as  those  used  in  the 
NTSC  video  standard. 

Other  examples  of  multichannel  2-D  signals  can  be  found  in  linear  arrays  when 
the  sensor  measurements  are  inherently  separated  along  certain  channels  (for  exam¬ 
ple  in-phase  and  quadrature,  principal  and  orthogonal  returns  of  a  radar,  and  in  the 
measurement  of  dual  polarized  radar  cross  section  as  a  function  of  space  and  time  for 
extended  targets.  In  these  examples  one  can  easily  imagine  multichannel  signals  that 
are  higher  than  two-  dimensional. 

Although  the  image  processing  community  has  long  dealt  with  multiple  correlated 
images,  [see  e.g.  Ref.  l],  this  community  has  for  the  most  part  not  approached  the 
analysis  from  a  signal  processing  point  of  view.  A  notable  exception  is  the  work  by 
Hunt  and  Kubler  [2]  who  seem  to  have  been  the  first  to  apply  signal-processing  methods 
to  multichannel  image  restoration.  The  signal  processing  community  on  the  other  hand 
has  dealt  with  the  analysis  of  multidimensional  signals  but  has  also  largely  ignored  the 
analysis  of  multiple  multidimensional  signals.  While  multichannel  2-D  signals  can  be 
regarded  as  a  special  case  of  3-D  signals,  multichannel  signals  have  their  own  special 
properties  that  makes  their  analysis  distinct  and  in  many  cases  more  tractible  than  the 
analysis  of  3-D  signals  in  general.  Therefore  we  have  chosen  to  concentrate  on  these 
special  signals  and  study  their  properties  explicitly. 

This  report  is  primarily  concerned  with  random  multichannel  2-D  signals  and  their 
applications.  These  signals  will  have  properties  sometimes  similar  to  multichannel 
1-D  signals  and  sometimes  more  like  2-D  signals.  We  will  be  concerned  here  mostly 
with  topics  related  to  filtering  and  estimation  of  these  signals;  and  in  particular  linear 
prediction  and  its  applications. 

The  remainder  of  this  report  is  organized  as  follows.  Chapter  II  deals  with  the  rep- 
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reservation  of  multichannel.  2-D  signals  and  their  statistical  characterization.  Chapter 
III  focuses  on  linear  prediction  and  autoregressive  (AR)  modeling.  This  forms  a  basis 
for  the  remainder  of  the  material.  Chapter  IV  discusses  spectrum  estimation.  Model- 
based  methods  for  estimating  the  entire  spectral  matrix  are  presented.  The  estimate 
includes  the  2-D  autospectrum  for  each  channel  and  the  magnitude  and  phase  of  the 
cross-  spectra.  Chapter  V  describes  applications  of  the  theory  to  image  processing. 
Chapter  IV  provides  a  summary  and  conclusion.  Areas  for  future  work  are  cited  there. 
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II.  Statistical  Characterization  of  Multichannel  2-D  Random  Processes 
2.1  Representation  of  Multichannel  2-D  Signals  and  Systems 

2.1.1  Signal  Domain  Representation 

Discrete  multichannel  2-D  signals  (or  sequences)  will  be  represented  by  a  vec¬ 
tor  quantity  x(n!,n2)  as  illustrated  in  Fig.  2.1.  Each  component  xK(n1,n2)  of  the 
vector  represents  a  2-D  signal  existing  in  one  of  the  channels.  Note  that  the  vector- 
valued  nature  of  the  function  relates  to  the  fact  that  the  signal  is  multichannel  while 
the  vector  nature  of  the  argument  n  =  (nj,n2)  relates  to  the  fact  that  the  signal  is 
multidimensional* 

The  processing  of  these  signals  is  analogous  to  that  for  other  discrete  signals.  We 
will  be  primarily  concerned  with  linear  shift-invariant  (LSI)  operations  that  can  be 
represented  by  2-D  vector  difference  equations.  In  particular  an  M-Channel  2-D  signal 
x(ni ,  n2 )  is  transformed  to  another  M-  channel  2-D  signal  by  an  LSI  system  represented 
by  the  difference  equation 

y(ni,n2)  =  —  AXli.y(n1-i1,n2-i2)+  Btlll  x(n1  -  lx ,  n2  -  l2 ) 

(U.ijjea 

(«i.»3)*(0,0) 

(2.1) 

where  the  Aili2  and  Blltl  are  M  x  M  matrix  coefficients  and  a  and  (5  define  their 
(fixed)  regions  of  support.  The  system  will  be  recursively  computable  if  there  exists  an 
ordering  for  processing  of  the  points  such  that  values  of  the  signal  y  needed  to  compute 
the  signal’s  current  value  at  the  point  (n1?n2)  are  always  available.  Whether  such  an 
ordering  exists  depends  on  the  shape  of  the  output  region  a. 

The  output  for  a  multichannel  2-D  system  can  be  equivalently  represented  by  the 
2-D  convolution  operation  in  either  of  the  forms 

CO  ,oo 

y(ni,n2)—  _  4)  (2.2a) 

f  -  CO  OO 

CO  ,  OO 

y(ni,n2)=  ^  H  (nx  -  ,  n2  -  £2)x(£i ,  t2)  (2.26) 

£  1  £  3  —  —  CO  ,  —  OO 

where  H (•,  •)  is  a  matrix  function  representing  the  2-D  multichannel  impulse  response. 
(Each  term  if(£l5£ 2)  in  Eq.  (2.2)  is  an  M  x  M  matrix.)  A  2-D  multichannel  system 


*  Although  it  is  tempting  to  use  vector  notation  for  the  arguments,  this  tends  to 
hide  the  essential  2-D  nature  of  the  signals.  Therefore  we  use  the  longer  but  more 
explicit  notation  involving  the  two  indices  nx  and  n2 . 
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Channel m 


x1(n),n2) 
X^,  ng) 

X3(nl  ,n2) 


Figure  2.1  Discrete  Multichannel  2-D  Signal. 


-  4- 


will  be  called  a  Finite  Impulse  Response  (FIR)  system  if  the  support  of  is  finite 

and  an  Infinite  Impulse  Response  (HR)  system  otherwise.  Clearly  for  an  FIR  system 
all  of  the  Aili7  in  Eq.  (2.1)  are  zero  and  H{tl,t2)  =  Blll2  for  (£i,£2)  €  /?. 

For  purposes  of  this  report,  we  will  be  dealing  with  statistical  properties  of  random 
signals.  Thus  the  mean  of  a  signal  x(n!,n2)  is  a  vector  quantity  defined  by 

m*  («i ,  n2 )  =  £[x(n! ,  n2 )]  (2.3) 


and  the  correlation  and  covariance  are  M  x  M  matrix  functions  defined  (respectively) 
by 

R(nl,n2  ;  m1,m2)  =  E  [x(n1,n2)xT  (m1,m2)]  (2.4a) 


C(n1,n2;m1,m2)  =  E  (x(n!,n2)  -  mjtinnj))  (x(m1,m2)  -  mI(m1,m2))J 


(2.46) 


When  the  random  process  representing  the  signal  is  homogenous  or  stationary  the  mean 
is  constant  and  the  correlation  and  covariance  are  functions  only  of  the  vector  distance 
between  the  points  (nl5n2)  and  (m1,m2).  In  this  case  Eqs.  (2.3)  and  (2.4)  become 


and 


m*  =  E  [x(nx ,  n2)] 


(2.5) 


R(ki ,  k2 )  —  E 
C(kl2k2)  =  E 


[x(n1,n2)xr(n1  -klyn2  -  k2)] 
(x(nj ,  n2)  -  mx )  (x(nx  -  k1 ,  n2 


-  k2)  -  m,)T 


(2.6a) 

(2.66) 


It  can  be  seen  from  their  definition  that  the  correlation  and  covariance  functions  have 
the  following  symmetry  properties 


R{kuk2)  =  RT{-k1,-k2)  (2.7  a) 

C{k1,k2)  =  CT{-k1,-k2)  (2.76) 


When  a  stationary  random  signal  is  transformed  by  a  linear  shift-  invariant  system, 
the  statistical  characteristics  of  the  output  can  be  expressed  in  terms  of  those  of  the 
-  input  by  using  Eqs.  (2.1)  and  (2.2)  and  taking  expectations.  Thus,  from  (2.1)  the 
mean  satisfies 

my  =  -  X  ^..>4™!/+  X] 

( *  x  i  *’  3 )  €  a  ( l  i  » * 3  )  £  P 

(*1  .*5  )t£  (0,0) 


or 


-  i 


T.  ^  •  1  ■  J 


E 


.  (  *  l  .  *  2  )  G  a 


mT 


(2.8) 
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where 


(2.9) 


The  correlation  function  satisfies  the  pair  of  equations 

Ry(k  i,k2)  =  —  Y,  Ail{lRy(k1 —i1,k2 —i2)+  Y,  Bll(2  Rxy(k1  —  ,  k2 —i2) 

(  *  1  1 1  2  )  €  Of  (^1^2)6/? 

( i  1 ,  *  2  )  7*  ( 0 , 0 ) 

(2.10a) 

Ryx{k  1  5  ^2)  “  ”  ^  ^  ^%i  i7  Ryx  (^1  ?  ^2  “*^2  )  “t*  ^  ^  R  £\ ,£3  (^1  “^1  5  ^2  ”^2  ) 

(*’  1  »**2  )  €  ot  (^i^2)€/? 

( *  1 » * 2  (0,0) 

(2.106) 

where  the  cross  correlation  is  defined  by 

Rxy{ki,k2)  =  Ryx{-k  1,-^2)  =  E  [(x(n1?n2)  -  ms)  (y(nx ,  n2)  -my)T]  (2.11) 
Alternate  relations,  in  terms  of  the  impulse  response  are,  from  (2.2) 


my 


H{i  i,4)  m3 


(2.12) 


and 

00 , 00  00  , 00 

Ry(k1,k2)=  Y  H{t1,l2)Rx(k1+p1-t1,k2+p2-t2)HT  (pi,p2) 

£  1 ,  £  2  —  —  00  ,  —  00  pi,p2  =  -oo,-oo 

(2.13) 

The  relation  between  the  correlation  and  the  cross  correlation  can  also  be  expressed  in 
terms  of  the  impulse  response  as 


OO  ,  OO 

Ryxi^i,^)-  Y  H(i1,t2)Rx(k1  -  £i,A;2  -  4)  (2.14a) 

(^1  »^2  )  =  —  00  -*  —  00 

and 

OO  ,00 

Ry(ki)k2)=  H{£\  5  ^2  )Rxy  (^1  ~  ^1  5  ^2  —  ^2  )  (2.146) 

£  l  :  £  2  —■  —  OO  ,  —  OO 

which  forms  a  connection  between  Eqs.  (2.10)  and  (2.13).  Relations  analogous  to  Eqs. 
(2.10),  (2.11),  (2.13),  and  (2.14)  clearly  apply  to  the  covariance  and  are  easily  derived. 

2.1.2  Frequency  Domain  Representation 

When  a  (deterministic)  signal  of  the  form 

x(nl5n2)  =  cz^z^  (2.15) 
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(2.16) 


is  applied  to  an  LSI  system,  the  output  (from  Eq.  2.2)  is  given  by 


y(nl5n2)  =  Hz  (z1 ,  z2)cz^  z” 5 


where 

co  ,  co 

Hz  (zx ,  z2)  =  ^  H{h,l2)z-tl  z;1’  (2.17) 

,  £2  =  “  OO  ,  —  OO 

is  the  z  transform  of  the  impulse  response  and  is  called  the  system  function.  If  one 
applies  Eqs.  (2.15)  and  (2.16)  to  Eq.  (2.1)  it  can  be  seen  that 


Hz  ( Zi ,  z2 ) 


E 


Aj  {  z, 
1 


.  (t'l  )€  a 


y1 

j 


E 


(2.18) 


The  frequency  response  of  the  system  is  the  system  function  evaluated  on  the  two  unit 
circles 

H,u{u1,oj2)=Hz  (2.19) 

Stability  of  the  2-D  multichannel  system  corresponds  to  convergence  of  Hz  on  the  unit 
circles.  Stability  thus  implies  that  the  frequency  response  exists  for  all  values  of  uq  and 
u>2  • 


Random  signals  are  characterized  in  the  frequency  domain  by  the  power  spectral 
density  matrix.  This  is  defined  as  the  2-D  Fourier  transform  of  the  correlation  function 


OO  .  CO 

■Muq,^)  =  Y,  Rx{kl,k2)e-j“'k'e~j“'-k'  (2.20) 

Jc  1  ,k  j  =  —  co,  —  00 


The  spectral  matrix  is  positive  semi-definite  and  Hermitian.  Its  off-diagonal  terms 
represent  cross  spectra  between  components  of  the  signal  and  need  not  be  real.  For  the 
case  of  a  two-channel  2-D  signal  the  matrix  can  be  written  as 


^(Wi,CJ2)  = 


Sn  (wi,w2) 

*^21  (W1  5  ^2  ) 


^1  2  (W1  1  W2  ) 

S22  (W1  5  ^2  ) 


(2.21) 


The  terms  Sxl  and  S22  are  real  and  nonnegative  and  represent  the  2-D  power  spectrum 
or  autospectrum  for  each  channel  while  the  terms  S12  and  S2l  represent  cross  spectra. 
Because  of  the  Hermitian  symmetry  the  cross  spectra  have  identical  magnitude  and 
opposite  phase.  The  normalized  quantity 


k 2(uq  ,u>2)  = 


l‘S'12  (wi  >  w2 )  |2 
■Si  !  (cUj  ,  u;2  )  .S22  (wi )  W2  ) 


(2.22) 


is  the  2-D  magnitude  squared  coherency  (MSC)  and  is  frequently  used  instead  of  |S'12  | 
to  represent  the  magnitude  of  the  2-D  cross  spectrum.  Further  treatment  of  the  MSC 
is  given  in  Appendix  A. 
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A  cross-spectral  matrix  Sxy  (uq ,  tu2)  can  also  be  defined  as  the  Fourier  transform 
of  the  cross  correlation  function 


(2.23) 


/c  x  ,/c  3  =  —  OO  ,  —  OO 


Its  components  represent  cross  spectra  between  components  of  the  signal  x  and  com¬ 
ponents  of  y.  This  matrix  is  neither  positive  definite  nor  Hermitian.  However  it  shares 
a  symmetry  property  with  Syx  in  that 


Sxy{oJi,oj2)  =  Sy*J  (wi,w2) 


(2.24) 


When  a  signal  is  transformed  by  an  LSI  system,  frequency  domain  expressions  for 
the  output  power  spectrum  and  the  cross  spectrum  between  input  and  output  can  be 
obtained.  These  are  derived  by  taking  Fourier  transforms  of  Eqs.  (2.13)  and  (2.14). 
The  output  power  spectral  matrix  has  the  form 

Sv  (w. ,  y2 )  =  H.v  (w. ,  ) S,  (Wl ,  uq )a:T(y1,u,)  (2.25) 

while  the  cross  spectral  matrices  relate  to  Sx  and  Sy  through  the  expression 

Syx  (wi )  (uq ,  c o2)Sx  (uq  ,  uq>)  (2.26a) 

Sy  (oq ,  u2 )  =  Hw  (u1,oj2)Sxy(u1,u2)  (2.266) 

and  the  symmetry  relation  (2.24). 

2.2  Vector  Representation  of  Multichannel  2-D  Signals 

2.2.1  Signal  Vector  Representation 

Frequently,  a  multichannel  signal  is  defined  over  some  finite  rectangular  region  of 
support.  That  is,  the  signal  is  considered  only  for  0  <  <  Ni ,  and  0  <  n2  <  N2. 

In  this  case  it  can  be  useful  to  represent  the  signal  as  an  NXN2M- dimensional  vector 
of  the  signal  values.  Two  such  representations  are  most  useful.  In  the  first,  the  vector 
valued  signals  at  points  ( nl,n2 )  are  lexicographically  organized  first  by  row  and  then 
by  column  into  the  larger  vector.  That  is,  the  multichannel  2-D  signal  is  represented 
by 


x  = 


*0 

*1 


(2.27a) 


where 


x(nl50) 
x(ni,  1) 


.x(n! ,  N2  -  1). 


(2.276) 
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This  representation  will  be  called  index  ordering. 


A  second  type  of  representation  organizes  signal  component  within  each  channel 
by  rows  and  columns  and  then  stacks  the  components.  This  representation  has  the 
form 


where 


x  = 


X  = 


-i 

^0 


— Ni  -  1  “J 


X1 


,M 


m  —  1,2, ...  ,m 


where 


K,0) 

(^1,1) 


ni  —  0|  1>  •  •  •  >  Ni  1 


(2.28a) 


(2.286) 


x"'  = 

■Im  (^1,-^2  -  1) 

This  representation  will  be  called  component  ordering. 

Index  and  component  orderings  are  related  by  permutation  transformations 


(2.28c) 


x'  =  Px 
X  =  Prx' 


(2.29a) 

(2.296) 


The  form  of  the  permutation  matrix  is  illustrated  in  Fig.  2.2  for  N1  =3 ,  N2  =  4,  and 

M  =  2. 


2.2.2  Statistical  Characterization  of  Signal  Vectors 

The  signal  vectors  previously  defined  can  be  characterized  by  mean  vectors  and 
correlation  or  covariance  matrices.  These  quantities  are  defined  for  index-ordered  vec¬ 
tors  by 


m  =  P[x]  (2.30a) 

R  =  F[xxr]  (2.306) 

K  =  P[(x  —  m)  (x  —  m)r  ]  (2.30c) 

and  for  component-ordered  vectors  by  similar  expressions  with  primed  variables.  From 
Eqs  (2.29)  and  (2.30)  it  can  be  seen  that  the  index-ordered  and  component-ordered 
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Figure  2.2  Permutation  matrix  for  Nj  =  3,  N2  =  4.  M  =  2  . 
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statistics  are  related  by 


m'  =  Pm  (2.31a) 

R'  =  PRPr  (2.316) 

K'  =  PKPt  (2.31c) 


Note  that  for  a  stationary  random  process  the  vector  m  will  consist  of  Nx  iV2  M- 
dimensional  subvectors  equal  to  the  signal  mean  (see  Eq.  2.5).  The  vector  m'  will 
consist  of  MiYx  iV2 -dimensional  subvectors.  Each  subvector  has  elements  that  are  all 
identical  and  equal  to  the  mean  of  the  signal  in  the  corresponding  channel. 

The  correlation  and  covariance  matrices  for  a  stationary  2-D  random  process  have 
specific  structure  and  symmetry  at  their  various  levels  of  partitioning.  These  properties 
are  illustrated  in  Figs.  2.3  and  2.4.  An  example  showing  detailed  structure  of  these 
matrices  is  given  in  Appendix  B. 

2.2.3  Linear  Transformation  of  Signal  Vectors 

Linear  transformations  on  signal  vectors  can  be  represented  by  matrix  equation 

y  =  Ax  (2.32) 

We  will  consider  only  index-ordered  transformations  here  since  index  ordering  is  more 
extensively  used  in  this  report.  The  equivalent  matrix  for  component-ordered  transfor¬ 
mation  is  given  by 

A'  =  PAPt  (2.33) 

General  linear  transformations  such  as  this  may  operate  upon  the  rows,  columns,  or 
channel  dimensions  of  the  signal. 


Linear  transformations  corresponding  to  a  LSI  filtering  operations  are  represented 
by  a  matrix  of  the  form 


H  = 


H(0)  H(— l) . . .  Ht-iVi+l) 

H(l)  H(0)  . . .  H(— iW+2) 


_H(iVa-l)  H(iW-2) 


where 


HK)  = 


iL(ni  ,0) 

H [nl ,  l) 


1) 

ffK,0) 


lH(ni,N2  —  l)  H{n1,N2-2) 


H(0) 


H{nlt-N,+  1) 
H (n: ,  —  iV2  +  2) 


(2.34a) 


(2.346) 
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Signal  Vector 


Xo 

x(n,.0) 

*1 

x(n,,l) 

where  xn  = 

— r*i 

XVx-1 

x(n,.;Y,-l) 

Correlation  Matrix 


RjO)  R(l) 

R(-l)  R(0) 


R=£  = 


R(A'i-l) 

RfV-2) 

(BLOCK  TOEPLLTZ  WITH 
AjA/xAVW  BLOCKS) 


R(  —  jV,  —  1)  R(-Ar,-2)  R(0) 


where 


R(k.O)  R(k.l)  .  R[k..\\-1) 

R(k,~  1)  R(k,0)  .  R(k..\\- 2) 


R(*)  =  £  xnxnr_* 


(BLOCK  TOEPLITZ  WITH 
M x M  BLOCKS) 


R(k.-.\,-l)  R(k- A'j-2)  R(Jk.O) 


where 


R(kl.ki)  =  RT(-kl-k,)  =  Ex(nl,n,)xT(rh-ki.n,-k2)\  (HOT  TOEPLITZ ) 

Fig.  2.3  Correlation  matrix  for  index  ordering. 


-  12- 


Signal  Vector 


X1 

-y  "1 
±0 

im('ll-O) 

o 

X" 

XT 

where  xm  = 

where  x™  = 
n\ 

xv 

v  rn 

X.Y.-l 

—  1 

Correlation  Matrix 


R.i 

R 12 

R 

R2J 

Ro2 

R 

R  '  =  £  x  'x  T 


(.V, A'.xAVVj  BLOCKS ) 


where 


pi  i c  ol  it  pi  it 

/t0  A!  ■‘SVj  - 1 

pi  k  ol  k  Dl  ^ 

n_x  n0  ...  rcVi_o 


R,  t  =  £  x'x*r  = 


(BLOCK  TOEPLITZ 
WITH  .Vo >  .V;  BLOCKS) 


pi  it 


i  A 

—  iV. 


/? 


i  it 
0 


where 


/? 


i  it 
n, 


=  F  v'x‘f 
-c.  x„xr>_„i 


( TOEPLITZ ) 


Fig.  2.4  Correlation  matrix  for  component  ordering. 
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(The  ~  notation  has  special  meaning  that  will  be  explained  shortly.)  When  the  trans¬ 
formation  corresponds  to  an  FIR  filter  many  of  the  blocks  and  subblocks  in  Eq.  (2.34) 
will  be  zero;  the  specific  arrangement  will  depend  on  the  shape  of  the  support  region 
for  the  impulse  response. 

For  a  system  described  by  a  linear  difference  equation  (2.1),  one  can  write 


Ay  =  Bx  (2.35) 


where  the  matrices  A  and  B  are  defined  in  terms  of  the  coefficient  matrices  and 

Bj  analogously  to  Eqs.  (2.34).  The  output  y  can  then  be  expressed  directly  as 


y  = 


(2.36) 


This  form  of  the  analysis  implicitly  assumes  zero  initial  conditions  for  the  difference 
equation. 

Occasionally  it  is  necessary  to  convert  between  one  form  of  channel  representation 
for  the  2-D  multichannel  signals  and  another  without  performing  any  specific  filtering 
or  signal  processing.  This  would  occur  for  example  if  the  signal  x(n.!  ,  n.2)  represented 
a  color  image  with  red,  green,  and  blue  components  and  we  wanted  to  convert  the 
signal  to  NTSC  form.  If  the  new  signals  are  defined  in  terms  of  the  old  signals  by  the 
transformation 

x(nltn2)  =  Tx(ni,n2)  (2.37) 

then  the  index-ordered  signal  vector  would  be  transformed  according  to 


x  =  Tx  (2.38) 

where  T  is  a  NlN2M  by  NXN2M  block  diagonal  matrix.  This  matrix  would  contain 
N !  N2  nonzero  blocks  all  equal  to  the  signal  transformation  T. 

2.2.4  Reversal  Notation 

In  the  analysis  of  2-D  signals  using  vector  representation,  it  is  sometimes  necessary 
to  deal  with  vectors  and  matrices  whose  components  are  ordered  corresponding  to 
decreasing  values  of  the  signal  index  rather  than  increasing  values  as  in  Eq.  (2.27). 
While  a  reordering  of  the  signal  values  can  be  represented  by  a  permutation  matrix 
operation,  it  is  much  clearer  to  develop  explicit  notation.  The  notation  to  be  used 
will  apply  only  to  index-ordered  vectors  and  matrices.  Component-ordered  quantities 
are  less  frequently  used  and  do  not  lend  themselves  readily  to  the  reorderings  that  are 
considered  here. 
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The  symbol  ~  above  a  vector  will  correspond  to  a  vector  whose  first-level  partitions 
have  been  organized  according  to  decreasing  value  of  the  index. 


*Wl-i 

Si/x-2 


*0 


(2.39) 


where  the  xn  are  defined  by  Eq.  (2.27b).  We  call  this  a  first  index  reversal  of  the 
vector.  The  symbol  ~  below  a  vector  will  be  used  to  denote  a  vector  whose  second 
level  partitions  have  been  organized  according  to  decreasing  index  values. 


(2.40) 


where 

"x^ ;  —  l)  - 

x(ni ,  No  —  2) 

x  = 

—  n’  1 

-  x(<i,,0) 


(2.41) 


We  call  this  a  second  index  reversal  of  the  vector  and  is  occasionally  useful.  Most 
frequently  it  will  be  required  to  use  the  doubly  reversed  vector  which  is  denoted  by  a 
double~over  the  vector: 


(2.42) 


& 


and  denotes  a  vector  where  both  the  main  blocks  and  the  second  level  blocks  have  been 
reversed. 


Since  transformation  of  a  reversed  vector  usually  implies  generation  of  a  newr  vector 
with  corresponding  order,  it  is  most  reasonable  to  define  reversals  of  a  matrix  as  a 
reordering  about  both  its  rows  and  columns.  Thus  a  matrix  A  represents  a  matrix  A 
whose  first  level  blpcks  or  partitions  have  been  reversed  along  both  rows  and  columns. 

Matrices  £  and  A  are  defined  similarly  with  reversals  about  rows  and  columns  of  the 
inner  blocks.  Note  also  that  reversal  about  both  rows  and  columns  of  a  matrix  is 
equivalent  to  transposition  about  both  the  main  and  reverse  diagonals.  This  preserves 
certain  symmetry  such  as  block  Toeplitz  structure  when  it  exists. 


If  reversal  notation  is  used,  the  LSI  transformation  of  Eq.  (2.2)  can  be  represented 


as 


y  =  Hx 


(2.43) 
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where  the  matrix  H  is  now  defined  by 


with 


H(nO  = 


H(0) 

H(l) 

...  H(iV1- 

H  = 

H(— 1) 

H(0) 

...  H(Ar1  — 

.H(— iVi  +  1)  H(— IVi  +  2) 

H(0) 

■  H(ny,  0) 

H{nlt  1) 

...  H{nltN2 

1) 

...  H  (riy ,  N2 

-H(rii,N2  —  1) 

H{ny,N2  -2) 

H (ni , 

(2.44a) 


(2.446) 


Reversal  notation  will  also  prove  to  be  useful  in  considering  signal  models  in  Chapters 
III  and  IV. 

2.3  Separable  Multichannel  2-D  Signals  and  Transformations 

2.3.1  Separable  Signals 

The  direct  product  of  an  L  x  Q  matrix  B  and  an  N  x  P  matrix  A  is  an  NL  x  PQ 
matrix 

B  cl  1 1  B  cl  1 2  ...  Bay  p 
Ba2l  Ba22  ...  Ba2p 


B  <8>  A  — 


(2-45) 


,  Ba^  i  Ba^j  2  ...  Ba ^  p  _ 

In  general  B  <g)  A  is  not  equal  to  A  0  B.  A  number  of  other  properties  are  given  in 
Table  2.1  and  in  Appendix  C. 


Table  2.1  -  Properties  of  Direct  Product 


(B  <S»  A)  (D  ®  C) 

=BD  <g>  AC 

{ B®A)T 

=BT  0  AT 

(B®A)~l 

=B~ 1  0  A~  1 

Consider  for  the  moment  a  single  channel  2-D  random  signal  that  has  the  separable 

form 

x(ni,n2)  =  xil(n1)  •  xB  (n2)  (2.46) 

Then  if  xA  and  xB  are  the  vector  representations  of  the  signals  (n^)  and  xB  (n2),  a 
vector  representation  of  the  2-D  signal  is  given  by 

X  =  xB  ®xA  (2-47) 
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Now  if  (nj )  and  xB  ( n2 )  are  independent  random  processes  with  correlation  matrices 
Ka  and  RB  the  correlation  matrix  for  the  2-D  signal  vector  is 


E  [xxT]  =  E  (xB  0  xA )  (xB  0  xA )' 

=  E  [(xbXb)]  ®  E  [(xAx^)] 


(2.48) 


—  RB  0  R>i 

A  similar  direct  product  relation  holds  for  the  mean  vector  and  the  covariance  matrix. 


Multichannel  2-D  random  processes  can  be  separable  along  the  index  (nx,n2) 
directions,  between  channels,  or  both.  Various  cases  of  separable  signals  and  their  cor¬ 
relation  matrices  are  given  in  Table  2.2  assuming  index  ordering.  The  specific  structure 
of  the  matrices  can  be  envisioned  by  comparisons  with  the  matrices  in  Fig.  2.3. 


Table  2.2  -  Separable  Forms  for  Multichannel  2-D  Signals 


Form  of  Signal 

Description 

Form  of  Correlation 

x(n.j ,  n2)  =  xA  (nx)  -xB  (n2) 

Product  of  a  1-D 
single-channel 
signal  and  a  1-D 
multichannel  signal 

x(nj,n2)  =  x(n1,n2)  •  c 

Product  of  a  2-D 
single  channel  signal 
and  a  M-dimensional 
random  vector 

R^  =  Rc  0  Rx 

x(n.j ,  n2)  =  xA  (nj  •  xD  (n2)  -c 

Product  of  two  1-  D 
single  channel  signals 
and  a  random  vector 

R^  =  Rc  0  Rb  0  R^ 

Separable  structure  for  the  signals  when  it  exists  leads  to  simplifications  in  the 
analysis.  For  example  2-D  problems  can  be  reduced  to  a  pair  of  1-D  problems.  The 
implications  of  separability  for  linear  prediction  are  described  in  Chapter  III. 

2.3.2  Separable  transformations 

If  a  linear  transformation  for  a  multichannel  2-D  signal  is  completely  separable  it 
can  be  represented  as  the  direct  product  of  three  matrices 

T  =  W  0  V  0  U  (2.49) 
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where  U  is  a  Ni  x  Ni  matrix  representing  transformation  along  the  nx  direction,  V  is  an 
N2  x  N2  matrix  representing  transformation  along  the  n2  direction  and  W  is  an  M  x  M 
matrix  representing  transformation  between  channels.  If  the  signal  is  correspondingly 
separable  as 

x  =  c®xB  ®x^  (2.50) 

then  the  processing  of  the  signal  is  greatly  simplified  since 

Tx  =  ( Wc )  ®  (FxB  )  ®  (UxA  )  (2.51) 

However  separable  transformations  are  advantageous  even  if  the  signal  is  not  separable. 
Common  examples  of  separable  transformations  are  the  DFT,  Hadamard  and  Walsh 
transforms,  and  the  singular  value  decomposition. 

A  convenient  relation  exists  between  index-ordered  and  component-ordered  sepa¬ 
rable  linear  transformations.  In  particular,  if  the  index-ordered  transformation  is  given 
by  (2.49)  then  the  component-ordered  transformation  is  given  by 

T  =  V  ®U  ®W  (2.52) 


A  somewhat  more  concise  representation  for  separable  transformations  exists.  This 
involves  first  representing  the  multichannel  2-D  signal  in  yet  another  form,  as  a  NXM  x 
N2  matrix 


r*i  i 


lxM] 


(2.53) 


where  Xm  is  a  x  N2  submatrix  whose  elements  are  the  signal  xm  [nl ,  n2)  in  channel 
m.  It  can  be  verified  that  with  this  representation,  a  separable  transformation  of  the 
form 

y  =  (W  ®  F  ®  17)  x  (2.54) 

leads  to  a  matrix  representation  Y  with  blocks 


M 

Ym  =J2^meUXeVT 
£=  1 


m  =  1,  2, . . . ,  M 


(2.55) 


where  wml  are  the  elements  of  W .  This  can  be  written  in  a  single  matrix  equation  as 


Y  =  (/  ®  W)  ■  UXVT 


(2.56) 
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III.  Linear  Prediction  for  Multichannel  2-D  Random  Processes 

3.1  Linear  Prediction 


3.1.1  Problem  Statement 

Let  x(n1,n2)  represent  a  zero-mean  stationary  multichannel  2-D  random  signal.  The 
linear  prediction  problem  is  concerned  with  forming  a  linear  estimate  of  the  signal  at  the 
point  (nl5n2)  from  other  values  of  the  signal  in  a  region  a.  Specifically  we  form  the 
estimate 

x(nl5n2)  =  -  ^  A?  i^{n1  —  h,n2  — 12)  (3.1) 

)€<* 

(«i,»j)*(0,0) 

and  the  prediction  error  is  defined  as 

e(n1,n2)  =  x(nl5n2)  -  x(nl5n2)  (3.2) 

The  matrix  coefficients  Aj  ;  are  chosen  to  minimize  the  mean-  squared  error  E  ' e ( n. j. .  )  | 3  . 

The  error  is  generated  from  the  data  by  the  FIR  prediction  error  filter.  From  Eqs. 
(3.1)  and  (3.2)  the  difference  equation  for  the  filter  is 

e(nl5n2)=  A J1<ix(n1  -  tL,n2  -  t2)  (3.3) 

(  *1  i  *  7  )  6  Of 

with  .Too  =  I ■  Since  the  filter  is  FIR,  the  terms  Aj^  .  also  represent  the  impulse  response 
and  a  is  the  impulse  response  region  of  support  (see  Eq.  2.2a).  The  correlation  function 
of  the  error  evaluated  at  lag  (0,0)  is  the  M  x  M  prediction  error  covariance  matrix  E,. 
That  is 

=  E  [£(n1,n2)er  (nl5n2)]  =  Re( 0,0)  (3.4) 

The  mean-squared  error  £  is  given  by 

£  =  E  [|c(n1,n2)|2]  =  E  [eT  (nx ,  n2)e(n1 ,  n2)]  =  trEe  (3.5) 

These  quantities  are  important  for  the  analysis  that  follows. 

3.1.2  Normal  Equations 

The  prediction  error  filter  coefficients  and  the  prediction  error  covariance  are  formed 
by  solving  a  set  of  Normal  equations.  These  linear  equations  follow  directly  from  the 
orthogonality  principle  [4]  which  states  that  the  error  must  be  orthogonal  to  the  signal 
values  used  in  prediction. 
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(3.6a) 


E  [x(nx  —  i'x ,  n2  -  t2)er  (^1 ,  n2)]  =  [O] 

(j'i,z2)€a:,  (*i,t2)  ^  (0,0) 

It  further  follows  from  Eqs.  (3.3),  (3.4)  and  (3.6a) 

E  [x(n1,n2)er(n1,n2)]  =  E  [e(nx ,  n2)er  (nx ,  n2)]  =  E£  (3.66) 

These  last  two  equations  can  thus  be  written  more  concisely  as 

E  [x(nx  —  fi,n2  —  f2)eT  (nl5n2)]  =  Ee6(i1)<5(i2) 

(i'i ,  f 2 )  G  Q!  (3.  / ) 


3. 1.2.1  Rectangular  Support 


The  Normal  equations  will  be  developed  by  using  an  index-ordered  vector  representa¬ 
tion  for  the  signal.  At  first  let  a  be  the  rectangular  region  shown  in  Fig.  3.1.  The  size  of  a  is 
Pl  xP2  points  and  Lx  and  L2  can  be  any  values  in  the  range  —Pi  <  Lx  <0,  — P2  <  L2  <  0. 
Then  an  index-ordered  representation  for  the  signal  data  that  occurs  in  Eq.  (3.3)  is  needed. 

To  make  the  analysis  clear,  consider  the  case  where  Ly  =  L2  =  0.  This  is  the  quarter 
plane  or  first  quadrant  predictor.  The  signal  data  needed  in  Eq.  (3.3)  is  in  the  range 
rix  —  Pi  +  1  to  n1  and  n2  —  P2  +  1  to  n2  (see  Fig.  3.2).  Form  a  vector  xni„2  as  follows 


wn  i  n  2 


-ni-P  i  +  l 

— n  i  —  P  i  +  2 


■^n  i  —  P  i  4*  * 


x[nl  -  Pi  +  t,  n2  —  P2  +  1) 
x(ni  -  Pi  +  i,  n2  -  P2  +  2) 


x(ni  -  Pi  +  t ,  n2 ) 
Then  Eq.  (3.3)  can  be  written  in  matrix  form  as 

where 


n2) 

=  Arx  „ 

— n  i  n 

r  a<°>  i 

— 

A(1) 

-4t0 

— 

-4a 

-■4>,f3-i  - 

(3.8a) 


(3.86) 


(3.9) 


(3.10a) 


(3.106) 


20 


‘2 


Figure  3.1  A  general  rectangular  region  of  support  ac 
for  the  prediction  error  filter. 
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I  a"  X(n-i 


i  ,=  0  V2 

l2~  0 


Figure  3.2  First  quadront  predictor. 
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The  Normal  equations  then  follow  directly.  From  Eqs.  (3.7)  through  (3.9)  we  have 


E 

xrllrt5€T(n1,n2) 

=  E 

or 


where 


=  E 


x„  „  x„  „ 

— n  i  n  i  — n  i  n  ; 


A  =  S 


(3.11) 


RA  =  S 


(3.12) 


rs(0>  i 


s  - 


o 


o 


(3.13a) 


(3.136) 


and  the  partitioning  in  Eqs.  (3.13)  corresponds  to  that  of  Eqs.  (3.8)  and  (3.10).  Eq.  (3.12) 
represents  the  Normal  equations.  The  results  for  a  first  quadrant  multichannel  predictor 
are  summarized  in  Fig.  3.3. 


Observe  that  the  correlation  matrix  that  occurs  when  the  Normal  equations  are  writ¬ 
ten  in  the  standard  form  (3.12)  is  the  reverse  matrix  R  and  not  R  (compare  Figs.  2.3 
and  3.3).  This  point  is  not  usually  appreciated  or  discussed  in  the  literature  on  linear 

prediction.  For  the  single  channel  case  R  and  R  are  both  the  same;  for  our  multichannel 
case  it  is  important  to  make  the  distinction. 


When  Lx  and  L2  of  Fig.  3.1  are  not  equal  to  zero  a  similar  analysis  again  leads  to 

Eq.  (3.12).  The  correlation  matrix  R  is  identical  to  the  one  for  the  previous  case  but  the 
matrix  of  filter  parameters  is  defined  by 


A  = 


A(Ll) 

A(Ll  +  Pl ~ 

4',l, 


A(,)  = 


-,d«,o 


.4,1) 


+  P  2  —  1  _ 


(3.14a) 


(3.146) 


23 


Prediction  Equation 


r  1  —  1/3—  i 

x(nLn2)  =  2  2  (i„t,H0.0) 

Ii  =  0  t2  =  0 


Normal  Equations 


R(0)  R(-l) 

R(l)  R(0) 

•  R(-P,-1) 

R(-P,-2) 

O  — 

.<<  ■ 

0  VL 

0 

R(P.-1)  R(P|-2)  • 

R(0) 

Aip.-" 

0 

where 


R{k,0) 

R(jt.-l) 

.  R(k,-P2-i) 

Se 

fl(M) 

R(k.O) 

.  R(k,-P2^  2) 

-Vl 

0 

R(k)=Rr(-k)= 

A<*>  = 

s(0|= 

R(k.P,-2)  . 

.  .  R{k. 0) 

Ak.pt-i 

0 

where 

( Atj . ^3)  =  T (-kl,-k2)- E\x(ni1n2)'KT ( A0  0=I  Ef=E  (x-x)(x-x)  r] 


Fig.  3.3  Equations  of  linear  prediction  for  first  quadrant  support. 
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The  matrix  S  is  similar  to  Eq.  (3.13)  but  the  non-zero  block  appears  in  a  position 
corresponding  to  the  location  of  the  coefficient  A0o  in  the  matrix  A. 

3. 1.2. 2  Special  Forms 

The  Normal  equations  for  rectangular  support  can  be  viewed  in  the  following  way. 

We  have  a  set  of  Pi  x  P2  points  as  shown  in  Fig.  3.4(a).  A  correlation  matrix  R  is  formed 
for  the  PXP2  points.  The  terms  of  the  correlation  function  appearing  in  this  matrix  are 
shown  in  Fig.  3.4(b).  The  Normal  equations  for  predicting  various  points  in  the  array  all 

involve  this  same  matrix  R.  They  differ  only  in  the  indexing  of  the  filter  coefficients  and 
in  the  position  of  the  terms  A0 o  and  in  the  arrays  A  and  S. 

The  filters  for  predicting  the  upper  right  and  upper  left  points  in  a  rectangular  region 
of  support  are  called  the  first  quadrant  and  second  quadrant  filters.  The  first  quadrant 
filter  was  discussed  in  detail  in  the  previous  subsection.  By  using  properties  of  the  reversal 
operation  the  second  quadrant  Normal  equations  can  be  put  in  the  form  of  Fig.  3.3.  When 
this  is  done  it  can  be  observed  that  the  equations  although  similar  have  their  second  level 
blocks  transposed.  Thus  the  filter  coefficients  for  the  first  and  second  quadrant  predictors 
are  not  the  same.  The  third  and  forth  quadrant  filters  for  predicting  the  bottom  left  and 
bottom  right  points  are  also  distinct.  Their  Normal  equations  differ  from  the  Normal 
equations  of  the  first  and  second  quadrant  filters  in  that  the  innermost  blocks  R(ki,k2) 
are  replaced  by  their  transposes.  This  is  a  difference  from  the  single  channel  2-D  case 
where  the  first  and  third  and  the  second  and  fourth  quadrant  filters  are  identical. 

A  very  important  form  of  predictor  is  the  nonsymmetric  half  plane  filter.  Although  this 
filter  clearly  has  non-rectangular  support  the  form  of  the  Normal  equations  can  be  derived 
from  those  already  considered.  The  linear  prediction  problem  is  depicted  in  Fig.  3.5.  One 
can  think  of  the  NSHP  support  region  as  a  rectangular  region  where  |L2|  points  (marked 
with  x’s  in  Fig.  3.5(a)  are  missing.  The  circled  point  is  the  one  predicted.  The  Normal 
equations  for  the  NSHP  can  be  obtained  by  starting  with  the  equations  for  rectangular 
support.  One  then  simply  drops  out  filter  coefficients  corresponding  to  the  “missing”  points 
and  eliminates  corresponding  rows  and  columns  of  the  correlation  matrix.  The  required 
terms  of  the  correlation  function  are  shown  in  Fig.  3.5(c).  The  detailed  form  of  the  Normal 
equations  is  given  in  Fig.  3.6. 

A  companion  NSHP  prediction  problem  is  shown  in  Fig.  3.5(b).  While  in  Fig.  3.5(a) 
data  is  processed  from  lower  left  to  upper  right,  in  Fig.  3.5(b)  data  is  processed  in  the 
opposite  direction.  We  refer  to  Fig.  3.5(a)  as  forward  prediction  and  Fig.  3.5(b)  as 
backward  prediction.  The  Normal  equations  for  the  backward  problem  can  be  obtained  in 
the  same  manner  discussed  for  the  forward  problem.  If  the  equations  are  permuted  and 
put  in  the  form  of  Fig.  3.6,  they  will  differ  only  in  that  the  inner  blocks  R(ki,k2)  appear 
transposed. 


25 


W-p,  -H 


>  n 


1 
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Figure  3.4  Linear  prediction  with  rectangular  support. 
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Note:  l_2  is  a  negative  number 

Figure  3.5  Linear  prediction  with  NSHP  support. 


Prediction  Equation 


h-P*-\  Pi- 1  Lt-Pt- 1 

x(n1n2)=  2  do^xlnj.n.-io)-  £  £  A*,*,3C( ni ~ *i*n2  — *2) 

1 2  =  0  1 1  ~  1  h~  L? 


Normal  Equations 


R  (o)  R'(-i) 

R'U)  R(°) 

.  R'(-P,  +  l) 

.  .  R(-P,+2) 

A'to) 

A(» 

S  '(0) 

0 

R'(P.-l)  RIP,- 2)  . 

.  RIO) 

A-"1'; 

0 

where 


R[k-L 2) 
R(k,-L2- 1) 

R(k.-Ln-l) 

R[k,-L2) 

R(k,~  L2-  P2-rl) 

.  R(k-L2-P2~ 2) 

Ak  ,L2 
Ak,L2+l 

R'(k)  =  R'T{-k)  = 

R(k:-L2-P2-l) 

R(k.-L2^P2- 2)  . 

R(k,  0) 

A(*)  = 

Ak,L2~P2-l 

and 


R  (o) = 

R(o,o) 

R(o.i) 

R(0.-1) 

R(0.0) 

.  .  R(0.-L2-P2-l) 

.  .  R(0,-L2-L2-P2^2) 

A'(°'  = 

O  O 

“  o 

S'(°)  = 

0 

R{Q.L2-P2-l) 

R(0.L2-  P2-2)  . 

R(0.0) 

ao,l2^p2-\ 

0 

and  where  all  other  quantities  are  defined  as  in  Fig.  3.3. 


Fig.  3.6  Equations  of  linear  prediction  for  NSHP  support  (forward  predictor). 
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3.2  Autoregressive  Models 


If  the  error  process  e(n.!  ,  n2)  produced  by  linear  prediction  were  white  then  one  could 
obtain  the  original  random  process  by  inverse  filtering.  Specifically,  by  replacing  e(n.l5n2) 
by  a  white  noise  w  and  solving  Eq.  (3.3)  for  x(n.!,n.2)  we  have 

x(n1,n2)  =  -  ^  Afi<3x(^i  -  i\,ri2  -  h)  +  w(nx,n2)  (3.15) 

( *  i .  *  3  )  €  a 

(*!  ,*3  )?S  (0,0) 

This  is  the  multichannel  2-D  autoregressive  (AR)  model.  The  white  noise  has  a  covariance 

T,w  =  E  [w(n1,n2)'wT  (nt,n2)]  (3.16) 

Fig.  3.7  illustrates  the  relation  between  linear  prediction  and  autoregressive  modeling. 
While  linear  prediction  uses  an  FIR  filter,  AR  modeling  employs'  a  recursive  or  HR  filter. 

Since  the  filter  coefficients  in  an  AR  model  satisfy  Normal  equations  identical  to  those 
for  linear  prediction,  there  is  a  temptation  to  equate  the  two  concepts.  The  distinction 
lies  in  the  fact  that  linear  prediction  does  not  usually  produce  white  noise.  When  one  uses 
an  AR  model  to  represent  an  arbitrary  random  process  and  sets  up  Normal  equations  to 
solve  for  the  filter  coefficients,  one  is  essentially  ignoring  the  distinction.  The  justification 
for  the  white  noise  AR  model  can  lie  only  in  the  belief  that  an  AR  model  of  some  order  is 
a  close  approximation  to  the  true  process.  For  a  model  with  NSHP  support  it  is  true  that 
the  AR  model  becomes  a  close  approximation  to  any  given  random  process  if  the  order 
becomes  sufficiently  large.  This  argument  follows  by  generalization  of  results  in  Refs.  5- 
and  6  to  the  multichannel  case. 


3.3  Linear  Prediction  for  Separable  Problems 

When  the  2-D  multichannel  random  process  is  separable,  various  simplifications  result 
in  the  linear  prediction  problem.  Because  the  correlation  matrix  is  separable  (see  Section 
2.3),  the  Normal  equations  (3.12)  for  linear  prediction  problems  with  rectangular  support 
reduce  to  a  set  of  lower-dimensional  Normal  equations.  The  results  for  a  predictor  with 
support  in  the  first  quadrant  are  summarized  here. 


When  the  correlation  within  channels  is  separable  from  the  correlation  between  chan¬ 
nels  the  covariance  matrix  R  in  (3.12)  can  be  represented  as  the  direct  product  of  a 

between-channel  covariance  matrix  Rc  and  a  within-channel  covariance  matrix  Rx .  In 
this  case  we  can  postulate  a  separable  form  for  the  matrices  A  and  S  and  show  that  they 
provide  a  solution  to  the  original  Normal  equations.  In  particular,  we  can  write 


(Rc  ®  R*  )  (I m  ®  a* ) 


(3.17) 


where  IM  is  the  M  x  M  identity  matrix,  a.x  is  an  N2 -dimensional  vector  whose  first 
component  is  1,  and  is  the  Ni  N2 -dimensional  vector  (l,  0, 0  . . .  ,  0)T .  Then  (3.17) 

implies  the  two  relations 


R*ax  =  a\tNx 


N, 


(3.18a) 
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(b)  AR  modeling 


Figure  3.7  Linear  prediction  and  A  R  modeling. 
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E(=o2zRc  (3.186) 

Since  the  solution  of  (3.12)  is  unique,  this  shows  that  it  can  be  found  from  (3.18).  The 
form  of  (3.17)  indicates  that  each  channel  has  the  same  filter  coefficients  ax  and  that  there 
are  no  between-channel  terms.  The  parameters  ax  and  a2  are  obtained  by  solving  the 
Normal  equations  (3.18a)  pertaining  to  a  single  channel  2-D  problem.  Then  from  (3.18b), 
the  prediction  error  covariance  for  the  2-D  multichannel  problem  is  a  scaled  version  of 
the  between-channel  covariance.  The  scale  factor  is  the  single  channel  prediction  error 
variance. 

Similar  results  can  be  obtained  if  the  correlation  matrix  is  separable  along  the  nl  and 
n2  directions.  In  this  case  the  correlation  matrix  can  be  represented  by  the  direct  product 
of  a  matrix  RB  which  relates  to  a  1-D  multichannel  process  along  the  n2  direction  and 
a  matrix  representing  a  1-D  single  channel  process  along  the  nl  direction  (see  Table 
2.2).  The  postulated  solution  is  of  the  form 


/  i  \ 

(R b  ®R*)(Ab  ®aA)  =  (  —  S,0)  1 

\a\  J 


& 


-.v. 


(3.19) 


which  implies 


RBAb  =  — S(0)  = 
oA 


Bp, 

0 

0 


(3.20a) 


R.4  (3.206) 

Equation  (3.20a)  is  a  set  of  Normal  equations  for  a  1-D  multichannel  problem  which  we 
solve  for  the  prediction  matrices  As  and  the  prediction  error  covariance  EPl .  Equation 
(3.20b)  represents  Normal  equations  for  a  1-D  single  channel  problem  that  we  solve  for  the 
filter  parameters  and  the  prediction  error  variance  aA .  The  solution  of  these  two  sets 
of  Normal  equations  then  allows  us  to  compute  the  filter  coefficients  as  As  ®  aA  and  the 
error  covariance  as 

=  aA  Ep,  (3.21) 

In  the  case  where  the  covariance  matrix  in  (3.12)  is  completely  separable,  along  rows 

and  columns  and  between  channels,  we  have  R  =  Rc  ®  RB  ®  R^  .  Then  we  are  led  to  the 
two  1-D  single  channel  subproblems 


R B  aJ3  —  °~B  L 


B 


and  the  2-D  multichannel  parameters  are  computed  from 

aA 


A  —  IM  ®  3b 

<?B  -n-C 


E,  =  o2  ol  R 


(3.22a) 

(3.226) 


(3.23a) 

(3.236) 
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The  decomposition  into  lower-order  linear  prediction  problems  discussed  here  would 
seem  to  be  restricted  to  the  case  where  the  multichannel  2-D  random  process  is  separable. 
In  actuality,  a  decomposition  of  the  multichannel  2-D  linear  prediction  problem  into  lower- 
order  problems  exits  in  general  where  there  is  no  separability.  This  is  discussed  in  the  next 
section.  A  difference  arises  however  in  that  the  lower-order  subproblem  for  the  general 
case  involves  an  expanded  multichannel  problem.  In  addition,  we  are  guaranteed  in  the 
separable  case,  that  if  the  correlation  matrix  in  (3.12)  has  doubly  block  Toeplitz  structure, 
the  correlation  matrices  of  each  of  the  subproblems  will  have  Toeplitz  or  block  Toeplitz 
structure.  In  the  general  case  the  Normal  equations  for  the  subproblems  do  not  all  involve 
correlation  matrices  with  Teoplitz  structure,  even  if  the  matrix  in  (3.12)  has  the  required 
doubly  block  Toeplitz  form. 

3.4  Relation  Between  Multidimensional  and  Multichannel 

Linear  Prediction 


There  is  a  close  relation  between  single  channel  2-D  linear  prediction  problems  and 
multichannel  1-D  linear  prediction  problems  (7j.  More  general  results  exist  that  show 
that  multidimensional  linear  prediction  problems  of  any  dimension  can  be  decomposed 
into  a  series  of  1-D  multichannel  and  single  channel  linear  prediction  problems  [8].  Some 
specific  results  will  be  discussed  here  that  show  how  2-D  multichannel  problems  relate  to 
higher  order  1-D  multichannel  problems.  These  results  will  be  used  in  the  next  section 
to  formulate  a  method  for  estimation  of  the  2-D  multichannel  model  parameters  without 
explicit  prior  estimation  of  the  correlation  function. 

Fig.  3.8  illustrates  a  multichannel  2-D  linear  prediction  problem  with  first  quadrant 
support  (a)  and  the  corresponding  multichannel  1-D  linear  prediction  problem  (b).  The 
Normal  equations  for  the  multichannel  2-D  problem  are  given  by  Eqs.  (3.12),  (3.13)  and 
(3.10).  The  structure  of  the  equations  is  further  detailed  in  Fig.  3.3. 

In  Fig.  3.8(b)  the  data  is  regarded  as  an  array  of  P2M  channels  of  1-D  signals  evolving 
along  the  nx  direction.  The  term  xn  used  in  the  index-ordered  representation  of  the  2-D 
data  corresponds  to  the  signals  in  this  array  of  channels.  The  1-D  P2 M-channel  linear 
prediction  problem  can  then  be  written  as 


Pi-  i 


(a<,))  2C.-1 


(3.24) 


»=  1 


where  xn  is  the  J°2 -^-dimensional  data  vector  and  xn  is  its  estimate.  The  coefficients  q:*'! 
have  the  form 


a(i) 
w0  0 

a(<) 

W01 

u0,Pj  — 

oc  ( 1 )  1= 

a**) 

U10 

( n ) 

U1  ,P2  - 

L  P 2  -1,0 

4*]- i.i  • 

..  a(i) 

P‘2  -  1,P2 
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Figure  3. 
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Multichannel  2-D  linear  prediction  and  related  multichannel  1-D  problem. 
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where  each  block  a.['}  is  of  size  M  x  M.  The  a(n)  are  found  by  solving  the  equation 


(3.26) 


where  EPl  is  the  P2M  X  P2M  prediction  error  covariance  and  where  the  correlation  matrix 
appearing  on  the  left  side  of  (3.26)  is  the  same  as  that  in  (3.12).  Now  post  mulitply  both 
sides  of  Eq.  (3.26)  by  the  term  Al0'  and  compare  the  result  to  Eq.  (3.12).  In  view  of  Eqs 
(3.10)  and  (3.13),  Eq.  (3.26)  will  be  identical  to  Eq.  (3.12)  if  we  require 

EP,A(0)=S(0)  (3.27) 

and 

A(i)=a(i)A(0)  ;  i  =  0, 1, . . . ,  —  1  (3.28) 

The  foregoing  equations  show  that  the  multichannel  2-D  Normal  equations  can  be 
solved  by  the  following  steps.  First  solve  the  P2 M-dimensional  1-D  multichannel  problem 
(3.26).  Since  this  problem  is  1-D,  the  multichannel  Levinson  recursion  can  be  employed 
to  find  the  coefficient  matrices  and  the  prediction  error  covariance  matrix  EPl .  Next 
solve  (3.27)  for  A*°E  this  is  also  a  set  of  1-D  Normal  equations  although  the  matrix  EPl 
is  not  in  general  Block  Teoplitz.  Finally  find  the  multichannel  2-D  coefficients  from  (3.28). 

Since  the  multichannel  2-D  coefficients  are  expressed  as  a  product  of  terms  in  (3.28), 
we  have  the  following  interpretation  of  the  multichannel  2-D  linear  prediction  problem. 
The  multichannel  2-D  prediction  error  can  be  computed  by  first  predicting  the  data  along 
the  direction  shown  in  Fig.  3.8(b).  This  is  a  1-D  linear  prediction  problem  and  results  in 
a  P2M-  dimensional  prediction  error  vector  corresponding  to  the  array  of  P2M  channels 
shown  in  the  figure.  This  prediction  error  is  itself  filtered  in  the  second  direction  to 
compute  the  M-dimensional  prediction  error  vector  for  the  2-D  problem.  The  Normal 
equations  solved  to  obtain  the  filter  coefficients  for  this  second  linear  prediction  problem 
are  represented  by  (3.27). 

This  view  of  multichannel  2-D  linear  prediction  has  one  further  interesting  aspect.  If 
the  multichannel  Levinson  recursion  is  used  to  solve  (3.26),  both  forward  and_  backward 
1-D  linear  prediction  parameters  are  generated  as  part  of  the  recursion.  By  arguments 
similar  to  those  presented  above,  the  backward  1-D  parameters  can  be  shown  to  relate 
to  the  coefficients  for  the  second  quadrant  2-D  filter.  Thus  an  algorithm  based  on  the 
multichannel  Levinson  recursion  can  be  used  to  solve  for  both  the  first  and  second  quadrant 
filter  coefficients  simultaneously.  This  will  be  shown  to  be  of  particular  use  for  the  spectrum 
estimation  problem  considered  in  Chapter  IV. 
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3.5  A  Direct  Method  for  AR  Parameter  Estimation 


The  connection  between  multidimensional  and  multichannel  linear  prediction  suggests 
a  method  for  estimating  the  AR  model  parameters.  This  method  will  be  called  a  direct 
method  since  it  does  not  require  prior  estimation  of  the  correlation  function. 

Our  method  capitalizes  on  the  existance  of  a  direct  method  for  solving  the  1-D  mul¬ 
tichannel  linear  prediction  problem.  We  will  refer  to  the  1-D  method  as  the  multichannel 
Burg  algorithm  since  it  is  based  on  ideas  originally  suggested  by  Burg  [9].  Details  of  the 
method  however  were  developed  separately  by  Nuttal  flOj  and  Strand  [  1 1] .  The  multi¬ 
channel  Burg  algorithm  is  described  in  Appendix  D. 

In  order  to  apply  the  multichannel  Burg  algorithm  to  2-D  data,  the  data  is  first 
partitioned  into  strips  along  the  n2  direction  as  shown  in  Fig.  3.9(a).  The  strips  of  width 
P2  for  an  Pi  x  P2  quadrant  filter  are  catenated  along  the  nj  direction  as  shown  in  Fig. 
3.9(b).  As  in  the  previous  section  this  data  is  considered  to  be  a  1-D  process  (in  the  n x 
direction)  with  P2M  channels.  The  multichannel  Burg  algorithm  is  then  used  to  estimate 
forward  and  backward  linear  prediction  parameters  (including  the  error  covariances).  This 
procedure  makes  no  explicit  estimate  for  the  correlation  matrix.  Discontinuities  where  the 
strips  were  catenated  together  are  ignored  since  they  represent  only  a  small  portion  of 
the  data.  The  error  covariance  matrices  that  are  computed  as  part  of  the  multichannel 
Burg  procedure  are  used  to  form  (3.27)  and  an  analogous  equation  for  the  backward  case. 
These  equations  are  solved  by  conventional  methods  and  the  multichannel  2-D  parameters 
are  computed  from  (3.28)  and  its  analog  for  the  backward  case.  This  completes  the  2-D 
estimation  procedure. 
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(b)  Multichannel  data  organization 


Figure  3.9  Sectioning  data  for  the  direct  method  of  AR  parameter  estimation 
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IV.  Multichannel  2-D  Spectrum  Analysis 
4.1  Spectrum  Estimation  Models 


Parametric  or  model-based  approaches  to  spectrum  estimation  are  based  on  a 
model  for  the  spectrum  involving  a  finite  and  relatively  small  number  of  unknown 
parameters.  These  parameters  can  be  estimated  from  the  data  and  used  in  a  formula 
that  gives  the  spectrum  of  the  model  in  terms  of  the  parameters. 

Autoregressive  modeling  is  one  form  of  spectrum  estimation  that  will  be  examined 
here.  A  model  of  the  form  of  Fig.  3.7(b)  is  used  to  represent  the  data.  Filter  parameters 
and  the  white  noise  covariance  are  estimated  by  the  techniques  described  in  Chapter 
III.  Since  the  input  spectrum  is  assumed  to  be  that  of  white  noise,  a  formula  based  on 
Eq.  (2.25)  can  be  used  to  compute  the  spectrum.  If 

H.(z,,z,)  =  Al.t,  (4-1) 

(£l  1^2  )  £ 


and  the  output  is  assumed  to  be  white  noise  with  constant  spectral  matrix  Sv^  ,  then 
the  spectrum  estimate  is  given  by 


S  (^! ,  u2 )  —  1  (^! ,  u2 )  Hu  T  (a?! ,  lo2  ) 

=  Hz  1  (zx  j  z2 )  Ew  Hj  (ax  1 ,  z2  1 ) 


zi  =  e 


(4.2) 


Note  that  this  results  in  an  estimation  of  the  entire  spectrum  matrix  i.e.,  all  of  the 
auto  spectra  and  cross  spectra  at  once  (see  Eq.  (2.2l)and  Appendix  A).  If  desired, 
normalized  quantities  such  as  the  magnitude  squared  coherence  can  be  computed  from 
the  elements  of  the  spectral  matrix  (see  Eq.  (2.22)  and  Appendix  A). 


Various  spectrum  estimates  can  be  obtained  by  assuming  different  support  for  the 
AR  model.  These  are  discussed  below. 


4.1.1  Non-svmmetric  Half  Plane  Models 

Non-symmetric  half  plane  models  are  a.  suitable  choice  for  the  spectrum  estima¬ 
tion  since  arbitrary  2-D  spectra  can  be  factored  into  sections  with  infinite  extent  non- 
symmetric  half  plane  support  [5].  Although  any  practical  algorithm  cannot  use  infinite 
or  even  very  large  support  regions,  nonsymmetric  half  plane  models  of  moderate  size 
have  proven  to  give  reasonable  results  for  spectrum  estimation.  An  interesting  fact  is 
that  it  is  unnecessary  to  consider  the  forward  and  backward  predictors  of  Fig.  3.5(a) 
and  (b)  spearately.  While  these  filters  have  different  parameters,  the  white  noise  covari¬ 
ance  is  related  such  that  formula  (4.2)  gives  theoretically  identical  spectrum  estimates. 
In  subsequent  sections  of  this  chapter  we  examine  results  of  NSHP  spectral  models  of 
various  orders.  The  definition  of  order  used  for  NSIIP  models  is  depicted  in  Fig.  4.1. 
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1  Definition  of  order  for  NSHP  models  . 
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4.1.2  Quadrant  Models 


AR  models  with  first  or  second  quadrant  support  are  of  interest  primarily  because 
of  their  convenience  of  computation  and  estimation  of  the  model  parameters.  The  third 
and  fourth  quadrant  filtes,  are  also  distinct  but  when  they  are  used  in  (4.2)  with  their 
corresponding  white  noise  covariance  they  give  estimates  identical  to  those  of  the  first 
and  second  quadrant  filters.  A  method  for  estimating  the  model  parameters  directly 
from  the  data  without  prior  explicit  estimation  of  the  covariance  function  was  described 
in  Section  3.5. 

It  will  be  seen  that  first  and  second  quadrant  AR  models  when  used  individu¬ 
ally  give  poor  estimates  of  the  spectrum  of  many  random  processes.  Features  of  the 
spectrum  such  as  peaks  tend  to  be  displaced  in  frequency  and  elongated  in  one  direc¬ 
tion.  However  a  certain  combination  of  the  models  has  proven  to  give  good  results 
for  estimation  of  spectra  in  2-D  single  channel  problems  (  Jackson  &  Chien  [12]  ].  We 
propose  here  a  generalization  of  this  method  to  2-D  multichannel  spectrum  estimation. 
A  combined  spectrum  estimate  is  computed  from 

S(wi,w2)  =  2  +  S,-/ (w:  ,w2))  1  (4.3) 

where  St  and  Su  are  the  spectra  corresponding  to  AR  models  with  first  and  second 
quadrant  support.  From  Eqs.  (4.3)  and  (4.2)  one  can  therefore  write 

S  [oji  ,u>2)  —  2  (^Hj  (uq  ,  u>2  )Dlv)  Hj  (uq  ,  w2)  +  Hn  (wi,w2  )Tiw1jj  Hj j  [u>i  ,u>2)^j  (4-4) 

In  the  special  case  that  the  noise  is  normalized  so  that  Ylw,  =  Ylw ,,  =  ^  Eq.  (4.4) 
has  the  simpler  form 

5(wi,w3)  =  2  {oj1,u2)Hi  (oq  ,w2)  +  (uq  ,  w2  )HU  (wx  ,  w2))  (4.5) 

This  form  is  analagous  to  the  form  proposed  by  Jackson  and  Chien  for  the  single  channel 
case. 

4.2  Resolution  and  Phase  Estimation  Experiments 

(Sinusoids  in  Noise  Background) 

Simulated  multichannel  2-D  random  signals  consisting  of  sinusoids  with  various 
relative  phases  and  frequencies  in  noise  backgrounds  were  generated.  This  section 
presents  estimated  spectra  of  these  random  fields. 

4.2.1  Comparison  of  NSHP  and  Quadrant  Modeling 

Here  we  discuss  the  estimation  of  a  sinusoid  in  noise  background  and  compare 
the  results  of  NSHP  modeling  to  quadrant  modeling.  Three  numerical  examples  are 
presented  which  illustrate  the  results. 
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Example  1 


This  example  is  concerned  with  the  analysis  of  spectra  for  a  two-channel 
2-D  single  sinusoid  signal  with  different  phase  in  additive  noise.  The  signals 
generated  in  channels  1  and  2  were 

xl  (nx ,  n2)  =  cos(n1w1  +  n2  w2)  +Wi  (%  ,  n2) 

x2  (n1 ,  n2)  =  cos(n1u>3  +  n2w4  +(/>)  +  w2  (n1 ,  n2 ) 

where  w1(n1 ,n2)  and  w2  (n,  ,  n2 )  are  zero-mean  independent  white  noise  sig¬ 
nals.  Spectrum  estimation  results  are  given  for  a  data  set  size  of  64  x  64, 
and  a  third  order  NSHP  filter.  Two  cases  are  considered  in  this  example.  In 
the  first  case  we  suppose  that,  the  two  channels  have  the  same  frequency,  but 
different  phase 


'^i  _  y,  =  -,u3  =  uq  =  -<p  =  1  radian j  , 

while  in  the  second  case,  we  assume  two  channels  with  different  frequency  and 
different  phase: 


=  u>2  =  —  ,u3  =  u>4  =  —  and  <f>  =  1  rad^  , 


Fig.  4.2  shows  the  results  for  the  components  of  the  2x2  spectrum  matrix  in  the 
first  case.  Only  the  cross  term  S12  (u>1  ,ui2)  is  shown  (magnitude  and  phase)  since  the 
term  S21(u1  ,u2)  is  theoretically  and  numerically  identical.  The  results  show  a  distinct 
peak  in  each  of  the  three  components  Su,S12,S22  corresponding  to  the  location  of 
the  sinusoid.  The  center  of  the  peak  is  accurately  located  near  (tt /2, tt/ 2).  Although 
the  phase  of  the  cross  spectrum  shows  various  artifacts  around  the  edge  of  the  region 
(where  the  magnitude  is  small  and  the  phase  is  that  of  the  noise)  the  phase  at  the 
location  of  the  sinusoid  is  accurately  estimated. 

Fig.  4.3  shows  5'(w1,w 2)  for  the  second  case  (different  frequency).  The  power 
spectrum  estimates  of  the  first  and  second  channels  *S'11(w1,w2)  as  S22(u1  ,u2)  have  a 
single-peak  at  the  location  of  the  sinusoid.  Although  the  cross  spectra  should  theoreti¬ 
cally  show  no  presense  of  sinusoids  some  small  amount  of  energy  is  detectable  at  those 
frequencies  in  the  cross  spectrum.  Similar  effects  have  been  observed  in  1-D  multichan¬ 
nel  spectrum  estimation  and  have  been  attributed  to  non  exact  pole-zero  cancellations 
in  the  estimate  for  the  cross  spectrum  [13]. 
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Figure  4.3  Estimate  of  spectra  for  sinusoids  of  different  frequency 
(NSHP  model). 
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Example  2: 

Two-channel  2-D  signals  with  two  sinusoids  in  additive  noise  are  considered 
in  this  example.  The  signals  in  the  channels  are  defined  by 

x i  (nj ,  ra2)  =  cos^u;,  —  n2 w2)  +  cos(n1  u>3  4-  n.2w4)  +  wl  [nl ,  n2) 
x2(nl,n2)  =  cos(n1w1  +  rc2u>2  +  4>i )  +  cos(niU;3  +  n2u>4  +  4>2) 

+  w2 (n.i ,  n2) 


The  results  are  presented  here  for  the  following  data:  u»,  (nx ,  n2 )  and  i v2[nl,n2) 
are  independent  zero-mean  white  noise  signals,  the  frequencies  and  phases  are 
u>i  =  u>2  =  =  cp 4  =  j,  and  <f>i  =  cf>2  =  1  radian  and  the  data  set  size  is 

64  x  64. 

The  power  spectrum  estimation  results  for  a  second  order  NSHP  model  are  given 
in  Fig.  4.4.  The  results  are  close  to  the  theoretical  results.  Su  (wi,o/2)  and  S22  (uq ,  u2 ) 
shows  that  the  two  sinusoids  are  easily  resolved.  The  estimated  amplitudes  are  unequal, 
but  this  characteristic  has  been  observed  in  even  1-  D  AR  spectrum  estimates.  The 
cross-spectrum  S12  (t<q  ,u>2 )  shows  the  sinusoids  resolved  and  the  phase  estimate  is  close 
to  the  true  phase  of  1  radian. 

Example  3: 

Three  sinuosids  in  each  channel  are  considered  in  this  example.  The  location 
of  the  sinusoids  is  shown  in  Fig.  4.5.  Observe  that  the  phase  of  the  sinusoids 
in  channel  2  differ  from  those  in  channel  1  by  1  radian.  The  sinusoids  are 
imbedded  in  white  noise  as  in  the  previous  examples.  Again  a  data  set  of  size 
64  x  64  was  used. 

Fig.  4.6  shows  the  results  of  spectrum  estimation  using  a  fourth  order  NSHP 
model.  The  results  show  good  estimation  of  the  position  of  the  sinusoids  in  the  au¬ 
tospectra  and  the  cross  spectra  and  there  is  almost  no  evidence  of  energy  from  sinusoids 
ujb  and  uiD  appearing  in  the  spectrum  of  the  opposite  channel  or  in  the  cross  spectrum. 
The  phase  of  the  cross  spectrum  in  the  region  where  the  sinusoids  are  located  is  nearly 
constant  with  a  correct  value  of  1  radian. 

Spectrum  estimates  for  Examples  1  and  2  were  also  computed  using  quadrant  plane 
models.  The  estimates  were  based  on  2  x  2  regions  of  support  for  the  first  and  second 
quadrant  filters.  Fig.  4.7(a)  shows  the  component  5X1  of  the  spectral  matrix.  Results 
for  S22  are  similar.  The  use  of  either  the  first  of  second  quadrant  filter  alone  results  in  a 
spreading  of  the  peak  in  one  direction.  The  combined  estimate  (Eq.  (4.3))  gives  a  more 
accurate  result  similar  to  that  of  the  NSHP  model.  Fig.  4.7(b)  shows  the  magnitude 
and  phase  estimates  for  the  cross  spectrum  Si2.  A  similar  spreading  phenomenon  is 
observed  in  the  magnitude  estimate  but  the  estimation  of  phase  is  correct  for  both  the 
individual  and  the  combined  spectrum  estimates. 
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Channel  1  tones 

Figure  4.5  Location  of  sinusoids  for  example  3. 
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Fig.  4.8  shows  the  spectrum  estimate  for  two  sinusoids  in  white  noise  corresponding 
to  Example  2.  In  this  case  the  phases  4>i  and  <+>2  were  taken  to  be  unequal  and  given 
values  of  1.5  and  0.5  radians  respectively.  The  results  of  this  quadrant  based  model 
can  be  qualitatively  compared  to  the  results  for  the  NSHP  model  of  Fig.  4.4.  The 
placement  of  the  sinusoid  along  a  diagonal  shows  some  characteristics  of  the  quadrant 
models.  The  first  quadrant  results  show  a  very  significant  spreading  of  the  peaks  along 
a  direction  orthogonal  to  the  line  connecting  their  centers.  This  is  observed  in  both 
the  autospectral  components  and  the  cross  spectra.  The  second  quadrant  estimates 
show  good  resolution  with  little  spreading  of  the  peak.  However  the  combined  estimate 
gives  the  best  results  with  sharp  peaks  and  with  magnitudes  more  nearly  equal  than 
those  observed  with  the  NSHP  model.  The  phase  in  Fig.  4.8(c)  is  slowly  varying  in  the 
region  of  the  sinusoids  for  all  of  the  estimates  with  a  correct  values  of  approximately 
1.5  and  0.5  radians  at  the  locations  of  the  sinusoids.  (Exact  values  produced  by  the 
combined  estimate  are  1.42  and  0.57  radians  respectively.) 

These  experiments  indicate  that  the  result  of  spectrum  estimation  using  a  single 
quadrant  model  is  not  generally  reliable  but  the  estimate  resulting  from  combining  the 
two  models  according  to  Eq.  (4.3)  is  quite  accurate. 

4.2.2  Model  order.  Data  Set  Size,  and  Signal-to-Noise 

Ratio  Experiments 

A  comprehensive  set  of  experiments  was  performed  to  determine  performance  of 
the  spectrum. estimation  procedures  as  a  function  of  model  order,  data  set  size,  and 
signal-to-noise  ratio.  The  results  are  briefly  summarized  here.  More  detailed  descrip¬ 
tion  of  the  results  will  be  reported  separately. 

The  results  of  the  model  order  and  data  set  size  experiments  showed  that  to  some 
extent  the  lack  of  resolution  resulting  from  a  small  data  set  size  could  be  compensated 
for  by  choosing  a  larger  model  order.  There  is  a  limit  to  this  trade-off  however  since 
a  larger  model  has  more  parameters  and  thus  should  require  more  data  to  estimate 
parameters  that  are  statistically  reliable.  The  experimental  observation  may  be  better 
restated  in  the  following  way.  When  the  data  set  is  large  a  smaller  order  model  can 
produce  results  that  are  comparable  to  a  large  model. 

The  case  of  two  sinusoids  in  noise  (Example  2)  was  repeated  for  smaller  size  data 
sets  using  a  NSHP  model.  For  a  second  order  filter  the  results  of  spectrum  estimation 
using  a  32  x  32  point  data  set  were  essentially  the  same  as  those  using  the  64  x  64 
point  data  set  reported  above.  For  smaller  data  sets  the  results  degraded  considerably. 
However  the  resolution  obtained  using  a  second  order  NSHP  filter  on  a  16  x  16  point 
data  set  was  similar  to  that  obtained  with  a  third  order  filter  on  an  '8  x  8  point  data  set. 
These  results  for  the  cross  spectrum  are  shown  in  Fig.  4.9.  Results  for  the  autospectra 
are  similar.  Note  in  Fig.  4.9  that  although  the  amplitude  of  the  estimate  varies  in  the 
three  cases  depicted,  the  phase  remains  essentially  constant. 
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Figure  4.7  Estimation  of  single  sinusoid  using  quadrant  models. 


1st  quadrant  2nd  quadrant 
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"Figure  *1.7  Estimation  of  single  sinusoid  using  quadrant  models.  (contM) 
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Fipure  4.8  Estimation  of  two  sinusoids  using  quadrant  models. 


FIRST  QUADRANT  SECOND  QUADRANT 
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Figure  ^.8  Estimation  of  two  sinusoids  us  .i  tip;  duadrant  models.  (eontM) 


FIRST  QUADRANT  COMB 
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Figure  ^.8  Estimation  of  two  sinusoids  using  quadrant  models.  (corit'd) 


Fip;ur.e  *1.9  KPfects  of  model,  order  and  data  set  size 


The  signal-to-noise  ratio  for  sinusoids  in  noise  was  defined  by 

SNR  =  10  log10  (4-6) 

where  A  is  the  amplitude  of  each  sinusoidal  component  and  o2  is  the  white  noise 
variance  (chosen  to  be  the  same  for  each  channel).  Fig.  4.10  shows  results  of  the  cross 
spectrum  estimation  for  a  64  x  64  data  set  with  closely  spaced  sinusoids  of  frequencies 
uq  =  u>2  —  7r/2  and  uq  =  u>2  =  the  model  used  for  these  experiments  was  second 
order.  At  a  SNR  of  12  dB  the  sinuosids  are  completely  resolved  with  sharp  peaks.  At 
a  SNR  of  3.5  dB  the  peaks  begin  to  merge  and  at  0  dB  the  peaks  become  a  single 
ridge  making  it  difficult  to  predict  that  there  are  two  sinusoids.  The  phase  obtained  by 
this  method  however  remains  essentially  constant  for  all  of  the  signal-to-noise  values. 
A  contour  plot  of  the  phase  is  shown  for  the  0  dB  case  again  illustrating  its  slowly 
varying  character  in  the  region  around  the  two  sinusoids. 

4.3  Estimation  of  More  General  Spectra 

This  section  contains  two  additional  examples  of  multichannel  2-D  spectrum  anal¬ 
ysis  involving  simulated  data.  These  cases  were  designed  to  test  the  accuracy  of  esti¬ 
mating  a  linear  phase  term  in  the  cross  spectrum  and  the  ability  to  estimate  parameters 
of  a  known  2-D  random  process  from  data.  The  latter  is  closely  related  to  the  2-D  lin¬ 
ear  system  identification  problem,  since  a  linear  system  can  be  identified  by  driving  it 
with  white  noise  and  then  estimating  the  cross  spectrum  between  input  and  output. 

4.3.1  Estimation  of  Delay 

For  this  experiment  two  channels  of  data  were  defined  by 
Zi  (n!,n2)  =  i ^(nj.na) 

x2  (n1,n2)  =  pXi  (nx  -  d2)  +  w2  [nl,n2) 

where  W 1  and  W2  are  two  independent  white  noise  processes.  Spectrum  estimates  were 
generated  and  the  slope  of  the  phase  of  the  cross  spectrum  was  measured  to  estimate 
the  delays  dx  and  d2.  Fig.  4.11  shows  the  cross  spectrum  estimate  that  was  obtained 
using  a  first  order  NSHP  model  for  white  noise  term  with  unit  variance  and  parameter 
values  of  (3  =  0.5,  d1  —  d2  —  1.  The  magnitude  of  the  cross  spectrum  is  constant 
and  the  phase  shows  a  linear  dependence  with  slope  corresponding  to  =  d2  =  1. 
Thus  the  experimental  result  agrees  precisely  with  the  theory. 

4.3.2  Estimation  of  Parameters  for  a  Linear  Model 

The  goal  of  this  experiment  was  to  estimate  the  parameters  of  a  known  linear 
process  from  records  of  data  generated  by  the  process.  The  following  multichannel  2-D 
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Figure  4.11  Linear  phase  estimation. 
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process  was  simulated. 
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where:  Wi(ni,n2)  (t  =  1,2)  are  two  independent  white  noise  signals.  This  process 
actually  has  support  only  in  the  first  quadrant  but  a  first  order  NSHP  model  was  used 
to  test  the  estimation  procedure.  For  this  problem  one  would  expect  the  two  matrix 
coefficients  Alt-  i  and  An  to  be  small  (near  zero). 

Table  4.1  shows  the  result  of  estimating  the  parameters  with  a  64  x  64  point 
data  set.  The  estimates  for  the  non-zero  parameters  of  the  model  are  close  to  the 
given  values,  and  all  of  the  remaining  parameters  but  one  are  at  least  one  order  of 
magnitude  smaller.  The  estimated  spectrum  components  for  this  process  are  shown  in 
Fig.  4.12.  Energy  is  spread  over  a  wide  range  of  frequencies  with  concentration  near 
higher  values  of  uq  and  lower  values  of  w2.  This  is  consistent  with  the  signs  of  the 
terms  in  the  defining  equations  for  the  process. 

Table  4.1  Estimated  parameters  for  a  1st  order  NSHP 
model  for  a  linear  process. 


^0,0 

1 
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A),l 

-0.56432 

-0.60426 

-0.027493 

-0.036813 

^-1,-  1 

0.014809 

0.00057926 

-0.14316 

-0.022275 

^1,0 

0.054647 

-0.012384 

0.58588 

-0.35247 

An 

0.017226 

-0.0-3484 

0.01443 

0.0055484 

4.4  Comparison  of  Direct  and  Indirect  Methods 

The  results  described  in  the  previous  sections  were  based  on  algorithms  that  first 
estimate  the  2-D  matrix  correlation  function  and  then  solve  Normal  equations  to  de¬ 
termine  the  AR  model  parameters.  These  methods  will  be  called  indirect  since  they 


-  57- 


0  2  0  4  06  0.8  1.0  0  0.2  0.4  0.6 


Figure  4.12 


Estimated  spectra  for  a  multichannel  2-D  linear  rrocess. 
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require  estimation  of  the  correlation  function  as  a  prerequisite  to  determining  the  model 
parameters.  A  direct  method  for  estimating  the  model  parameters,  i.e.,  a  method  that 
estimates  the  parameters  directly  from  the  data  was  described  in  section  3.5.  This 
method  is  very  well  suited  to  spectrum  estimation  using  combined  first  and  second 
quadrant  models  since  it  estimates  both  sets  of  filter  parameters  simultaneously.  Al¬ 
though  the  direct  method  has  not  been  exhaustively  tested  on  all  of  the  cases  reported 
earlier,  our  initial  results  indicate  that  its  estimates  are  at  least  comparable  to  and  in 
some  cases  better  than  those  of  the  indirect  methods. 

Figure  4.13  shows  the  results  of  estimating  a  single  sinusoid  in  white  noise  using 
combined  2x2  point  quadrant  filters.  These  results  are  the  same  as  the  combined 
results  shown  in  Fig.  4.7  but  are  depicted  in  a  slightly  different  format.  Fig.  4.14 
shows  the  spectrum  estimation  results  for  the  same  data  using  the  direct  method  of 
parameter  estimation.  The  direct  method  yields  a  sharper  peak  with  somewhat  lower 
sidelobes.  Estimates  of  the  one  radian  phase  shift  between  the  channels  are  1.05  rad. 
for  the  indirect  method  and  1.03  rad.  for  the  direct  method. 

Figure  4.15  shows  the  results  for  estimating  two  sinusoids  in  white  noise  using  com¬ 
bined  3x3  point  quadrant  filters  with  parameters  estimated  by  the  indirect  method. 
These  results  are  the  same  as  those  shown  earlier  in  Fig.  4.8.  Fig.  4.16  shows  the 
results  of  the  direct  method  applied  to  the  same  data.  In  this  case  the  results  appear 
to  be  nearly  equivalent.  A  few  minor  peaks  appear  in  both  cases.  For  the  indirect 
method  these  peaks  appear  closer  to  the  main  peak  while  for  the  direct  method  they 
appear  further  out.  Phase  estimates  at  the  peak  locations  are  1.47  and  0.55  radians 
for  the  indirect  method  and  1.48  and  0.50  radians  for  the  direct  method.  The  values 
produced  by  the  direct  method  are  slightly  closer  to  the  true  values  of  1.5  and  0.5 
radians. 
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Figure  4.13  Estimation  of  sinusoid  in  white  noise; 
indirect  parameter  estimation. 
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Estimation  of  sinusoid  in  white  noise; 
indirect  parameter  estimation,  (cont'd) 
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Figure  4.14  Estimation  of  sinusoid  in  white  noise; 
direct  parameter  estimation. 
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direct  parameter  estimation,  (cont'd) 
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Figure  b.15  Estimation  of  two  sinusoids  in  white  noise ; 
indirect  parameter  estimation. 
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Figure  4.15  Estimation  of  two  sinusoids  in  white  noise; 

indirect  parameter  estimation.  (cont'a) 


-  65- 


0.2  0.4  0.6  0.8  1.0  AWPLfTUDE 


1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 

♦ 

❖ 

J _ I _ I _ I _ I _ t  —  i-  i  t 

0  0.2  0.4  0.6  0.8  1.0 


u1 


<N 

3 


Figure  4.16 


(a)  s 

11 

Estiiration  of  two  sinusoids  in  white 
direct  parameter  estimation. 


(b) 


noise ; 


S 


22 


-  66- 


AWPLfTUDe 


(c)  s 

12 

Figure  4.16  Estimation  of  two  sinusoids  in  white  noise; 

direct  parameter  estimation.  (cont'd) 
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V.  Image  Processing  Applications 

This  chapter  deals  briefly  with  application  of  the  theory  of  multichannel  2-D  signals 
to  image  processing.  Three  applications  are  discussed  namely,  segmentation  of  color 
images,  image  coding,  and  image  spectrum  analysis. 

5.1  Image  Segmentation 

This  section  describes  an  algorithm  for  image  segmentation  that  uses  a  multichan¬ 
nel  2-D  model  to  represent  the  images  involved.  This  type  of  model  works  best  for 
images  involving  multiple  textures,  such  as  aerial  photographs  of  the  ground  in  rural 
areas. 

The  application  here  will  be  to  color  or  “false  color”  images  (the  latter  are  derived 
from  images  recorded  on  color  infrared  film).  For  these  cases  the  multichannel  2-D 
signal  has  three  components  (M=3)  corresponding  to  the  red,  green,  and  blue  video 
signals  (see  Fig.  5.1).  The  same  methods  could  be  applied  to  image  data  from  a 
multispectral  scanning  satellite  that  records  data  on  four,  seven,  or  even  more  IR 
channels.  In  this  case  the  number  of  components  M  would  be  equal  to  the  number  of 
channels  of  scanner  data  employed. 

It  is  assumed  that  regions  with  similar  texture  are  defined  by  boundaries  of  corre¬ 
sponding  pixels  in  each  of  the  three  color  image  components  and  that  within  each  region 
the  texture  is  described  by  a  multichannel  2-D  AR  model.  A  superimposed  Markov 
model  characterizes  the  occurance  of  regions  by  specifying  a  form  for  the  probability 
that  a  given  pixel  belongs  to  the  region  given  that  neighboring  pixels  belong  to  the 
same  or  other  regions.  The  combination  of  the  two  models  has  been  called  a  doubly 
stochastic  image  model  in  analogy  with  doubly  stochastic  models  occuring  in  the  study 
of  point  processes. 

The  doubly  stochastic  image  model  has  been  used  in  conjunction  with  monochrome 
images  and  single-channel  AR  models.  Since  many  of  the  details  of  the  algorithms  are 
described  in  our  earlier  work  [14]  for  monochrome  images  our  description  of  the  method 
here  will  be  brief. 


5.1.1  Segmentation  Algorithm 


The  image  to  be  segmented  is  assumed  composed  of  several  homogenous  regions 
of  texture.  The  red,  green,  and  blue  intensities  are  represented  by  the  multichannel 
2-D  signal* 


F(ni,n2)  = 


Fi{ni,n2) 

F2  (nL  ,n2) 
F3  (nx ,  n2 ) 


(5.1) 


* 


We  use  the  notation  F  to  represent  an  image  to  adhere  to  previous  convention. 
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g.  5.1  Representation  of  a  color  image  as  a  3-channel  2-D  random  process. 
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and  within  a  particular  region  this  signal  is  generated  by  the  AR  model 

F' (nx ,  n2)  =  Av.F'K  -  *1  ,n2  ~  *2)  +  W(nl5n2)  (5.2a) 

*1*2 

F  (nl9n2)  =  Ff(nl9n2)  +  p  (5.26) 

where  /x  is  a  constant  mean  vector.  The  white  noise  W(nl5n2)  is  described  by  a  3  X 
3  spatially  invariant  covariance  matrix 

=  E  [W(nx  ,n2)Wr  (rci,rc2)]  (5.3) 

which  in  general  is  not  diagonal.  A  separate  model  of  this  type  is  formed  for  each  of 
the  image  regions  in  which  there  is  a  different  texture. 

The  essence  of  the  segmentation  algorithm  is  as  follows.  Using  the  models  (5.2) 
and  a  Gaussian  white  noise  source,  we  can  form  a  probability  density  function  for  the 
set  of  all  pixels  in  the  color  image  conditioned  on  the  regions.  Assume  that  there  are 
Q  regions  R1 ,  £23  •  •  - ,  Rq  •  Since  the  texture  is  independent  from  region  to  region,  this 
density  function  can  be  written  as 


Q 

p(F|^1,^a,...,^0)  =  JIp(F|^,.)  (5.4) 

t=  1 

where  p  (F|£j)  represents  the  joint  density  for  the  pixels  contained  in  region  £,  .  Now 
assume  that  there  is  a  way  to  represent  the  probability  of  occurence  of  a  particular  set 
of  regions  in  the  image  Pr  [Z1,  R2, . . . ,  ZQ\.  Then  Bayes’s  rule  states  that 


Pr{R1,R2,...,RQ\F] 


p(F|£i£2)---i^q)  -Pr  [#! ,  £2 , . . .  Rq\ 

P(F) 


Thus  a  maximum  a  posteriori  (MAP)  estimate  for  the  regions  can  be  obtained  if  the 
Ri  are  chosen  to  maximize 


P  { F|£i5£2,...  Rq)  Pr  [Ri  y  R2)  - -.Rq] 


(5.5) 


Since  the  form  of  the  probability  density  function  for  the  white  noise  is  known, 
the  form  of  the  density  function  for  pixels  within  a  given  region  can  be  obtained  by 
transformation.  The  resulting  image  probability  density  function  is  most  conveniently 
represented  in  terms  of  the  errors  of  linear  prediction  which  can  be  written  as 

F(n.j ,  n2)  =  F  [nx ,  n2 )  —  ^  ^  t  i2  F  (n^  —  ,  n2  —  z2)  (5-6) 

*i.*i 
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where  the  assumed  constant  mean  level  is  subtracted  out  prior  to  applying  (5.6).  By 
using  the  linear  model,  the  density  function  for  pixels  within  a  region  can  be  repre¬ 
sented  as 

p(F|£,)=  Yl  P»k,  (E(ni,n2))  (5.7) 

where  pWk  is  the  probability  density  function  for  the  white  noise  source  of  type  k{ 
within  region  £,  . 

Now,  suppose  that  an  image  has  many  regions,  but  that  each  region  contains  only 
one  or  another  of  two  texture  types.  Then  choosing  a  set  of  regions  for  the  image 
is  equivalent  to  labeling  the  pixels  with  labels  0  and  1.  If  an  appropriate  statistical 
model  for  the  labels  of  the  pixels  exists  then  the  prior  probability  Pr  ,  £2,  •  •  •  £q] 
can  be  computed.  This  is  the  other  ingredient  needed  to  form  the  MAP  estimate.  An 
appropriate  statistical  model  is  one  that  represents  a  particular  class  of  2-D  Markov 
processes  [14,  15].  This  model  permits  calculations  of  the  desired  prior  probabilities. 

Combination  of  the  linear  filtering  model  with  the  Markov  model  results  in  an 
algorithm  to  obtain  a  MAP  region  estimate.  Since  the  resulting  equations  are  of  high 
dimension  and  nonlinear,  a  (possibly  suboptimal)  solution  is  obtained  by  iterating  the 
conditions 

[E<°»(n1,n2)]T[EW0]-1E<°>(n1,n2) 

0 

+^|SWo  I  -  2lnPr  [0|S(ni,n,)]  > 

1 

E(1)(ni,n2)]T  E(1)(nx,n2) 

+/n|Eu,I|-2/nPr  [l|S{rei,„3)]  (5.8) 

which  are  derived  from  (5.4)  through  (5.7)  and  the  Gaussian  form  of  the  density  func¬ 
tion.  The  terms  El°*(n!,n2)  and  E(1l(n1?n2)  are  the  error  terms  at  pixel  (nl5n2) 
computed  using  the  filters  of  type  0  and  1  respectively  and  Ea,.  is  the  correspond¬ 
ing  white  noise  covariance.  The  terms  Pr  [o| S(„ x >ri 3 ) ]  and  Pr  [l|5(raiin2)]  are  Markov 
transition  probabilities  representing  the  probability  that  the  pixel  at  location  (nl5n2) 
has  label  0  or  1  given  that  a  set  of  other  pixels  S(„i  n3)  in  the  neighborhood  has  a 
prespecified  set  of  labels.  The  inequality  provides  a  rule  for  assigning  pixel  labels  indi¬ 
vidually  based  on  labels  obtained  at  a  previous  iteration;  the  two  expressions  in  (5.8) 
are  evaluated  and  the  pixel  is  given  a  label  corresponding  to  the  sense  of  the  inequality. 

5.1.2  Experimental  Results 

Results  of  applying  the  segmentation  to  color  images  are  given  here.  Fig.  5.2(a) 
shows  a  128  X  128-pixel  color  image  representing  an  aerial  photograph  of  the  ground 
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(c)  green  component 


(d)  blue  component 


Figure  5.2  Color  image  of  trees  and  field. 
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in  a  rural  area.*  The  green  textured  area  is  a  grove  of  trees  and  the  predominantly 
yellow  area  is  a  field  with  tall  grass.  Figs.  5.2(b)  and  (c)  and  (d)  show  the  red,  green, 
and  blue  components  of  the  image.  (The  blue  component  is  very  dark  and  may  not 
reproduce  well.) 

Figure  5.3  shows  the  results  of  segmentation  of  this  image.  The  algorithm  used 
3-channel  2x2  first  quadrant  models  for  each  of  the  two  terrain  types.  The  Markov 
model  used  a  neighborhood  size  of  5x5.  Note  that  attempts  to  segment  this  image 
based  on  color  would  have  failed  because  the  fields  contain  large  linear  patterns  of 
green.  The  linear  filtering  models  however  match  the  colored  texture  in  the  image  and 
result  in  a  segmentation,  that  except  for  a  few  isolated  pixels,  is  very  accurate. 

Figure  5.4(a)  shows  another  color  image  of  trees  and  a  field.  This  image  is  some¬ 
what  more  difficult  to  segment  because  the  areas  of  green  in  the  field  do  not  appear  as 
a  linear  pattern  but  are  mixed  in  more  homogeneously.  Figure  5.4(b)  shows  the  results 
of  segmentation  using  image  models  similar  to  those  used  in  the  other  example.  Again, 
with  the  exception  of  some  small  areas  the  segmentation  is  quite  accurate. 

5.2  Image  Coding 

The  models  developed  in  this  report  can  be  applied  to  the  coding  of  color  or  other 
multichannel  images  for  purposes  of  storage  or  transmission.  In  this  application,  a 
model  is  formed  for  the  image,  the  parameters  of  the  model  are  stored  or  transmitted, 
and  the  model  is  used  to  reproduce  the  image.  Two  general  schemes  appear  to  be 
practical. 

The  first  scheme  is  applicable  when  portions  of  an  image  each  have  a  homogenous 
texture  and  one  is  concerned  about  reproducing  only  the  general  character  of  the  texture 
and  not  the  specific  gray  levels  at  each  pixel.  In  this  case  a  white  noise-driven  linear 
model  can  serve  to  characterize  the  image  sufficiently  so  that  only  the  parameters  of 
the  model  need  to  be  retained.  This  is  analogous  to  the  analysis-synthesis  method  for 
speech  encoding  [16] .  In  a  practical  application  of  the  method,  the  image  would  first 
be  segmented  and  the  boundaries  of  the  regions  would  be  coded  and  retained.  Within 
each  region,  a  linear  model  could  be  formed  for  the  image  and  the  parameters  (filter 
coefficients  and  noise  covariance)  would  be  estimated  and  retained.  The  image  would 
be  “decoded”  by  driving  the  filter  for  each  region  with  white  noise  having  the  specified 
covariance  and  thus  reproducing  the  image  texture  within  that  region. 

The  second  scheme  is  quite  general  and  would  apply  to  any  image.  It  is  analogous 
to  differential  pulse  code  modulation  (DPCM)  in  speech.  A  linear  predictive  filter 
is  derived  and  applied  to  the  image.  The  coefficients  for  the  filter  and  the  actual 
error  residuals  are  retained  and  coded.  Since  there  is  usually  great  redundancy  in  a 
multichannel  image  the  error  residuals  will  tend  to  have  low  dynamic  range  and  can 

*  Data  courtesy  of  Rome  Air  Development  Center,  Griffiss  AFB,  N.  Y. 
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Figure  5.3 


Segmentation  of  the  tree-field  color  image 
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(b)  segmentation 

Figure  5.4  Color  image  of  another  area  of 

trees  and  field  and  segmentation. 
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therefore  be  coded  with  a  fewer  number  of  bits  than  the  original  image.  The  image  is 
reproduced  by  inverse  filtering  using  the  decoded  error  residual  as  input  to  the  filter. 

The  method  just  described  is  most  effective  when  the  linear  predictive  filter  is 
well  matched  to  the  image  being  coded.  Since  one  cannot  expect  that  large  general 
images  will  be  well-represented  by  a  linear  filtering  model,  it  is  advisable  to  divide  the 
image  into  small  rectangular  blocks  and  develop  separate  linear  predictive  filters  for 
each  block.  Since  the  image  within  each  block  is  usually  more  homogenous,  one  can 
expect  that  lower  error  residuals  will  be  produced,  resulting  in  more  efficient  coding. 

DPCM  methods  have  been  developed  and  used  in  the  coding  of  single  images  and 
image  sequences.  Their  use  for  multichannel  images  using  the  linear  predictive  models 
discussed  in  this  report  would  seem  to  be  particularly  effective. 

5.3  Image  Spectrum  Analysis 

Spectra  for  portions  of  the  color  image  of  Fig.  5.2  were  computed  using  the 
methods  of  Chapter  IV.  Only  the  red  and  green  components  were  analyzed;  these 
were  designated  as  channel  1  and  and  channel  2  respectively. 

Two  portions  of  the  image  were  selected:  a  64  x  64-pixel  section  (one  quarter  of  the 
image)  in  the  lower  left  representing  the  field  and  a  64  x  64-pixel  section  in  the  upper 
right  representing  the  trees.  The  signal  index  nx  was  taken  as  the  horizontal  direction 
and  n2  was  taken  as  the  vertical  direction.  The  images  showed  distinct  differences  in 
their  spectra  which  could  be  used  as  a  basis  for  discrimination. 

Fig.  5.5  shows  the  2-D  spectral  components  of  the  fields.  The  red,  green,  and 
magnitude  of  the  cross  spectrum  appear  to  be  very  similar.  Power  is  distributed  in  the 
middle  and  higher  frequency  ranges  in  the  uq  direction  and  around  zero  frequency  in 
the  cj2  direction.  This  power  distribution  can  be  attributed  to  the  somewhat  evenly 
spaced  lines  of  green  (perhaps  some  kind  of  plants)  appearing  in  rows  in  the  field.  The 
phase  of  the  cross  spectrum  seems  to  oscillate  between  approximately  10°  and  -40°  in 
a  more-or-less  regular  fashion. 

Fig.  5.6  shows  the  estimated  spectral  component  for  the  trees.  Again,  the  red, 
green,  and  magnitude  of  the  cross  spectrum  are  very  similar.  The  power  tends  to 
concentrate  in  a  lower  region  of  the  cjx  scale  than  the  power  for  the  fields  with  a  peak 
occuring  around  7t/6  corresponding  to  the  spatial  placement  of  the  trees.  The  phase 
of  the  cross  spectrum  again  varies  between  +10°  and  -40°  but  with  somewhat  milder 
oscillations. 

The  concentration  of  power  around  u;2  =  0  in  the  tree  spectra  at  first  seemed 
strange.  One  would  expect  that  the  tree  image  would  show  similar  intensity  variations 
in  the  horizontal  and  vertical  directions  and  so  the  spectra  should  be  more-or-less 
similar  in  the  cjx  and  u>2  directions.  On  closer  examination  it  was  found  that  the 
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Fig.  5.5  Spectra  for  red  and  green  components  of  the  fields  image. 
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Fie  5.G  Spectra  for  red  and  trreen  comnoneiits  of  the  trees  irnaere. 


image  intensity  variations  for  this  data  are  not  the  same  in  both  horizontal  and  vertical 
directions.  Fig.  5.7  shows  typical  slices  of  the  image  in  the  nx  and  n2  directions  (we 
have  examined  several).  The  intensity  exhibits  rapid  changes  in  the  nj  (horizontal) 
direction  but  shows  a  very  slow  change  in  the  n2  (vertical)  direction.  This  accounts  for 
the  energy  in  the  spectra  being  concentrated  around  cu3  =  0. 

We  can  only  guess  a  reason  for  this  lack  of  symmetry  in  the  tree  image  data.  A 
plausible  explanation  is  that  a  slight  blurring  of  the  image  occurred  due  to  the  motion 
of  the  aircraft  during  the  time  at  which  the  photograph  was  taken.  This  would  have 
resulted  in  a  lack  of  higher  frequency  detail  in  the  direction  of  motion. 
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Fig.  5.7  Spatial  variation  of  intensity  in  the  trees  image. 


VI.  Conclusions 


This  report  developed  analysis  methods  for  multichannel  2-D  random  signals.  Al¬ 
though  these  signals  could  be  considered  as  a  special  case  of  three-dimensional  signals, 
the  multichannel  nature  of  their  origin  provides  them  with  special  properties  that  are 
not  shared  by  multidimensional  signals  in  general.  A  good  example  of  this  is  the 
concept  of  phase  of  the  cross  spectrum.  This  idea  seems  natural  when  considering 
a  multichannel  2-D  signal  but  would  have  no  counterpart  in  a  3-D  analysis.  While 
there  has  been  a  large  body  of  research  for  two  and  three-dimensional  signals,  analysis 
methods  specifically  for  multichannel  2-D  signals  have  not  heretofore  been  given  much 
consideration. 

A  substantial  part  of  the  work  in  this  report  concentrated  on  2-D  linear  prediction 
and  autoregressive  modeling.  Concepts  from  the  analysis  of  both  single  channel  2-D 
signals  and  multichannel  1-D  signals  were  generalized  in  this  development. 

An  important  application  of  the  signal  modeling  arises  in  spectrum  analysis.  The 
methods  described  in  Chapter  IV  of  this  report  showed  howautoregressive  models  could 
be  used  to  estimate  the  entire  spectral  matrix  for  multiple  2-D  signals.  Estimates  are 
developed  simultaneously  for  the  autospectrum  of  each  signal  and  the  magnitude  and 
phase  of  the  cross  spectra  between  the  signals  in  each  channel.  Spectrum  estimation  for 
multichannel  2-D  signals  is  the  topic  of  a  separate  investigation  [17]  and  many  results 
in  addition  to  those  reported  here  have  already  been  developed. 

The  application  of  the  analysis  methods  to  image  processing  was  discussed  briefly. 
Methods  for  image  segmentation,  image  coding,  and  spectrum  analysis  were  discussed 
in  this  report.  Many  further  image  processing  applications  such  as  enhancement  filter¬ 
ing,  target  detection  [18],  and  others  seem  possible. 

While  the  report  represents  the  results  of  a  serious  effort  to  develop  and  apply 
the  theory  of  multiple  2-D  random  signals,  the  work  is  by  no  means  complete.  Several 
applications  to  image  processing  have  already  been  mentioned  and  should  be  developed. 
In  addition,  the  analysis  of  some  types  of  radar  and  array  data  seem  to  be  important 
applications  of  the  spectrum  analysis  and  modeling.  Special  sensors  such  as  the  laser 
radar  [19]  where  the  data  collected  represents  intensity,  range,  and  doppler  in  two- 
dimensions  offer  further  opportunities  for  research. 

By  concentrating  on  linear  prediction  and  AR  modeling,  we  have  looked  at  only 
a  single  (although  rich)  aspect  of  the  analysis  of  multichannel  2-D  signals.  Other 
more  general  methods  of  modeling  and  estimation  including  pole-zero  and  noncausal 
modeling  were  not  fully  developed.  Representation  of  the  signals  in  other  ‘coordinate’ 
systems  by  applying  transformations  between  channels  may  lead  to  data  reduction  or 
other  advantages  in  various  applications.  We  hope  not  only  to  continue  our  work  in 
this  area  but  also,  by  publication  of  this  report,  to  stimulate  interest  in  the  further 
analysis  of  multichannel  two-dimensional  signals. 
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Appendix  A 

Spectral  Representation  of  Multichannel  2-D  Random  Processes 


A  stationary  multichannel  2-D  random  process  can  be  characterized  in  the  fre¬ 
quency  domain  by  a  2-D  spectral  matrix.  The  spectral  matrix  is  defined  by  Eq.  (2.20) 
of  Chapter  II  and  has  the  specific  form 


Sx  (c^i  ,w2) 


Si  1  (W1  )  ^2  ) 

5’i2(wi,w2) 

_  SM  !  (c^i  ,  0J2  ) 

5m2(Wi,W2) 

$1M  (^1  5  ^2  ) 

Sm  M  (tUl  ,  w2) 


(A.  1) 


where  S<7- (u^  ,  w2 )  is  the  2-D  Fourier  transform  of  Ri3{kl,k2)  the  cross  correlation  be¬ 
tween  channels  i  and  j.  The  spectral  matrix  is  a  periodic  function  of  uq  and  u2  with 
period  27r  in  each  dimension  and  so  it  needs  only  to  be  considered  only  on  the  interval 
—  IT  <  (Jj1  <  7T  and  —  7T  <  U)2  <  7T. 


It  follows  from  Eqs.  (2.7a)  and  (2.11)  that 

S,(wi,w2)  =  Sf  (Wi,U>2) 

and  thus  the  spectral  matrix  is  Hermitian  symmetric.  In  addition,  since  the  compo¬ 
nents  Ra  [kl ,  k2)  are  positive  definite  functions,  the  diagonal  terms  Si{  (w!  ,u>2)  which 
represent  the  2-D  power  spectrum  of  each  channel  are  non-negative.  Cross  spectral 
terms  Si3  (o^  ,  oj2  )  for  i  ^  j  are  complex  in  general.  However,  since  Ri}  (Aq ,  k2 )  is  real 
the  magnitude  of  Si}-  is  an  even  function  and  the  phase  is  an  odd  function  of  cj1  and 
C^2  * 


It  is  sometimes  convenient  to  write  the  matrix  of  the  squared  magnitudes  of  the 
terms  of  as 


Sn  ...  0 


0  ...  5 


MM  J 


^  2  ^2 


2  *  *  *  ^1M  fSn 


«!<,  Kht  •••  1 


0  1 


?M  M  J 


where  the  middle  quantity  is  called  the  squared  coherency  matrix.  It’s  elements  are  the 
magnitude  squared  coherency  (MSC)  between  the  it/l  and  jt/l  channels  and  are  defined 
as 


«?i(w  i,w2) 


Su{ui,U2)S]3(u>  i.,U>2) 


(A.  2) 


It  follows  from  the  properties  discussed  above  that  the  matrix  is  real  and  symmetric, 
that  is,  (w! , u>2 )  =  ,w2).  The  MSC  measures  the  correlation  between  the  the 

signals  in  channels  i  and  j  at  frequencies  ui3  and  u>2.  It  is  insensitive  to  spatial  shifts 
between  the  signals  in  the  two  channels.  These  effects  are  manifested  as  phase  in  the 
cross  spectrum. 
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Appendix  B 

Examples  of  Correlation  Matrices 
This  appendix  gives  examples  of  the  correlation  matrices 


R  —  E  [x  xT] 


and 


R '  =  E  x'x' 


for  index-ordered  and  component-ordered  representations  of  a  stationary  multichannel 
2-D  signal.  The  examples  are  given  for  Ni  =  3,N2  =  4,  and  M  =  2.  The  terms  r^e 
appearing  in  the  matrices  are  defined  by 


r'l-  =  E  [xi  (nx ,  n2 )xy  (nx  -  k,n2  -  t)} 


That  is,  r’lj’  is  the  element  in  the  ith  row  and  the  jth  column  of  the  matrix  correlation 
function  R(k,l)  (see  Eq.  2.6a). 
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Note:  r..  =  r.. 

'I  Ji 

The  matrix  is;  Block  Toeplitz  but  not  block  symmetric,  has  block  Toeplitz 
(but  not  block  symmetric)  blocks.  The  sub-blocks  are 
symmetric  but  not  Toeplitz  in  general. 

Figure  B.l  Index-ordered  correlation  matrix  =  3,  N2  =  4,  M  =  2* 
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Each  major  block  is  identical  except  for  subscripts  on  r.  The  overall  matrix 
is  neither  block  Toeplitz  nor  block  symmetric.  Major  blocks  are  block  Toeplitz 
(but  not  block  symmetric)  with  Toeplitz  blocks. 

Figure  B.  2  Component -ordered  correlation  matrix  N]  =  3,  N2=  4, 

M  =  2. 
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Appendix  C 

Direct  Product  of  Matrices 


The  direct  product  or  Kronecker  product  of  matrices  is  an  important  notational 
operation  in  dealing  with  multidimensional  signals.  This  appendix  summarizes  the 
most  important  properties  of  the  direct  product  without  proof.  For  a  more  thorough 
discussion  see  Ref.  [20]. 

Given  two  matrices 


A  = 


an  a12  •••  cl1p 


and 


aN 1  aN 2 


^ii  b i ; 


B  = 


aN  P 


bu 


(C.l) 


(02) 


_bLl  bL2  ...  bLQ  _ 

the  direct  product  B  ®  A  is  defined  as  the  partitioned  matrix 


B  ®  A  = 


Abi  i 

Abi2 

Ab^Q 

Abi,  x 

AbL  2 

. .  AbLQ 

(03) 


If  matrices  A,  C  and  B,  D  are  conformable  then  some  properties  are  as  follows. 


[B  ®  A)  {D  ®  C)  =  BD  ®  AC  (C.4) 

If  A  and  B  are  square  matrices  [N  =  P,  L  =  Q)  then 

(. B  ®  A)T  =  BT  ®At  (<7.5) 

tr(B  ®  A)  =  tr{B)tr{A)  (C-6) 

\B  ®  A\  =  \B\n \A\l  (C. 7) 

Further  if  A  and  B  are  nonsingular  then  the  direct  product  is  nonsingular  and 

{B  ®  A)'1  =  B'1  ®  A'1  '  (C. 8) 


Some  further  useful  properties  follow  from  the  preceeding  ones.  If  A,$  and  0,^? 
are  the  eigenvalue  and  eigenvector  matrices  for  B  and  A  respectively,  then  A  00,$®^ 
are  the  eigenvalue  and  eigenvector  matrices  for  B  ®  A.  That  is 

{B  ®  A)  ($  ®  <&)  =  ($  ®  (A  0  0)  (C.9) 
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To  state  this  another  way,  if  4>i  and  are  eigenvectors  of  B  and  A  corresponding 
to  eigenvalues  A,  and  then  4>i  0  if),-  is  an  eigenvector  of  B  0  A  corresponding  to 
eigenvalue  A That  is, 


(B  ®  A)  ■  {fa  0  Tpj)  -  A i9j(6i  0  Vv) 


(C.10) 


In  addition  if  A  and  B  are  symmetric  with  a  modified 


Cholesky  decomposition 


A  =  GDGt 


and 

B  =  G'D'G ,T 

where  G  and  G'  and  unit  lower  triangular  and  D  and  D'  are  diagonal,  then 


(B  0  .4)  =  (G'  ®  G)(D'  0  D)[G'  ®  G')T  (C.lil 


Some  further  fundamental  properties  and  cautions  should  be  observed.  It  is  clear 
from  the  definition  that  the  direct  product  is  associative,  i.e., 

(C  0  B)  0  A  =  C  0  (B  0  A) 

In  addition,  the  operation  is  distributive  over  addition: 

(B  +  C)®A  =  B®A  +  C®A  (CU3) 

A  0  (5  4"  C)  =  ^405  +  410^7  (C.14) 

However,  the  direct  product  is  not  commutative  (in  general  B  0  A  40  B)  and  it  is 
not  distributive  over  ordinary  matrix  multiplication. 


(c.  12) 


-  87- 


Appendix  D 

Multichannel  Levinson  Recursion  and  Burg  Algorithm 
Multichannel  Levinson  Recursion  (l-D) 

The  multichannel  Levinson  recursion  (also  called  the  Levinson-Wiggins-Robinson 
(LWR)  algorithm)  deals  with  a  linear  prediction  problem  of  the  form  (3.24)  and  solves 
the  Normal  equations  (3.26).  The  associated  backward  prediction  problem  is 


if  1-1 


x'  =  -  V  (a{i)'Yx 

■X-n-Ni  +  l  N  x  +  l  +  i 


>  =  1 


and  leads  to  the  Normal  equations'* 


R  ct'P1  =  R 


I 

p^i 

a(l)' 

= 

0 

_  CE1  "  ‘  -  1  > '  . 

0 

(3.24)' 


(3.26)' 


where  a*’)  has  a  form  analogous  to  (3.25).  The  recursion  is  specified  by  the  following 
set  of  equations. 


A,  =  [Rr  (l)  RT  (2)  ...  RT({)  ]  d^_!  (£>.la) 

a;  =Af=[R(l)  R(2)  ...  »(■)]«;„  (fl.lt) 

r<  =  (E;.,)"A,  (£>••  2a) 

r;  =E1-,1Af  (£>.26) 


a,  -i " 

0 

0 

A-i. 

0 

r, 

(D.3a) 

0 

r; 

{D.3b) 

e,  =  e,_,  (i-r;r,)  (dao) 

e;  =  b;_.  (i-r,r;.)  (£>.46) 

*  Note  that  ^  denotes  second  level  block  reversal. 
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The  terms  I\  and  T'  are  the  forward  and  backward  reflection  coefficient  or  partial  cor¬ 
relation  matrices.  The  term  A,  (A')  can  be  interpreted  as  a  cross  correlation  between 
the  forward  (backward)  prediction  error  and  the  one  unit  delayed,  backward  (forward) 
prediction  error.  Eqs.  (D.l)  -  (D.4)  are  applied  recursively  beginning  with  i  —  1  and 
ending  with  i  =  Px  —  1. 


Multichannel  Burg  Algorithm 


This  multichannel  form  of  the  Burg  algorithm  was  developed  by  Nuttal  f  10]  and 
Strand  [ll].  The  procedure  centers  around  estimating  the  terms  A,  and  A'  of  the 
multichannel  Levinson  recursion  directly  from  the  data.  Then  the  rest  of  the  equations 
(D.2)  -  (D.4)  can  be  applied  without  modification. 

Define  the  forward  and  backward  error  terms  for  the  1-D  multichannel  linear  pre¬ 
diction  problem  as 

§1*’  =x„  -x„  =  ajx  {D.5a) 

e'j  0  =  xn  _  ,•+  !-*!_<+!=  oc[T  x  (D.56) 

From  these  definitions  and  Eqs.  (D.3)  it  can  be  seen  that  the  errors  satisfy  the  recursion 


P(0  -  P(‘-i) 

Sn  —  in 


rfe'J-D 


*  —n 
T 


ell<)  =e'(i),  -r  ~(0 


— n  •+■  1 


t  %i  + 1 


(D.6a) 

(D.66) 


The  recursion  is  sustained  by  defining 

ri=  1 

n  =  1 


[D.l  a) 
( D-7b ) 
(D.7c) 


where  NP  is  the  number  of  points  in  the  data  record  minus  the  length  of  the  filter. 
The  matrix  A{  is  then  obtained  as  the  solution  of  the  bilinear  equation 


Cx 

C2 

C3 


+  a,c2  =  c3 

(£>.8) 

=  (e;)_1b 

(D.9a) 

=  (E  0“lD 

(D.9b) 

=  — 2G 

(D.9c) 
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where 


There  are  several  approaches  to  solving  (D.8).  One  approach  is  to  transform  A,  into 
an  equivalent  vector  representation  by  concatenating  rows  of  the  matrix  into  the  MNX- 
dimensional  vector  6.  The  matrix  C3  is  likewise  transformed  to  a  vector  c3 .  Equation 
(D.8)  can  then  be  written  as 

Q6  =  c3  (D.10) 

where 

Q  =  Cx  <g>  Iw  +  IM  <g>  C2  (D.ll) 

where  1^  is  the  M  X  M  identity  matrix.  Eq.  (D.10)  can  then  be  solved  for  6. 
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Principal  Symbols 


x(nl5n2) 

Multichannel  2-D  signal 

Zm  (nl5n2) 

Signal  component 

x,xn 

Signal  vectors  -  index  ordering 

5  A  ) 

Signal  vectors  -  component  ordering 

x,xm 

Signal  matrices 

AT  .  ,P.T. 

*1»J  5  *1  ‘2 

Coefficients  of  difference  equation,  filter  coefficients 

Impulse  response 

Hz  (2i  ,  22  ) 

System  function  (z  transform  of  impulse)  response 

Hw  (u^Wj) 

Frequency  response  (Fourier  transform  of 
impulse  response) 

H,  HK) 

Matrix  representation  of  impulse  response 

M 

number  of  channels 

m 

Mean  for  signal  x 

Rx  (A^  ,  k,2  )  ,  i?z  y  (^j  ,  k2 ) 

Correlation  functions 

Kx(k1,k2),Kxy(kl,k2) 

Covariance  functions 

Sx(u1,u2),Sxv(u1,w2) 

Power  density  spectra 

Kx  (Wi,W2) 

Coherence  function 

m 

Mean  vector  ,  index  ordering 

m' 

Mean  vector,  component  ordering 

R,  R(n) 

Correlation  matrix  index  ordering 

R',Rninj,Rnk'n' 

Correlation  matrix, 
component  ordering 

2, 

Multichannel  2-D  prediction  error  covariance  matrix 
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